Microsoft Word - 3-PETI#6580 21-31.docx Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 Free Vibration of Thick FGM Plates under TSDT and Thermal Environment Chih-Chiang Hong * Department of Mechanical Engineering, Hsiuping University of Science and Technology, Taichung, Taiwan Received 15 October 2020; received in revised form 30 November 2020; accepted 30 December 2020 DOI: https://doi.org/10.46604/peti.2021.6580 Abstract Three parameters of thermal environment, varied calculated shear correction, and third-order shear deformation theory (TSDT) of displacement are important in the frequency study. These three effects have been studied on the non-dimensional and dimensional frequencies of thick FGM plates. An additional c1 displacement term in nonlinear coefficient of TSDT is used to present the frequency of vibration into the simply homogeneous equation of thick FGM plates. The determinant of the coefficient matrix containing the c1 displacement term in dynamic differential equilibrium equations can be derived into the five degree polynomial free vibration equation. The non-dimensional and dimensional of natural frequency can be obtained. The effects of plate thickness, temperature of environment and power law index of FGM on the non-dimensional and dimensional frequency of FGM plates are investigated. Keywords: TSDT, FGM, frequency, nonlinear, homogeneous 1. Introduction Recently, many vibration papers devoted to the studies of mechanics of homogeneous, functionally graded material (FGM) porous beams/plates at micro, macro and nanoscales, diverse plates/beams theories, and linear/nonlinear thickness changes. In 2020, Żur et al. [1] presented the free vibration by using the nonlocal sinusoidal higher-order shear deformation theory (HSDT) for the magneto-electro-elastic FGM nanoplates. In 2020, Ouakad et al. [2] presented the nonlinear vibration by using the stress-driven nonlocal integral elasticity for the actuated hybrid nanotubes. In 2020, Ghobadi et al. [3] presented the nonlinear free vibration by using the Galerkin’s and perturbation methods for the FGM flexoelectric nano-plates. In 2020, Ghobadi et al. [4] presented the nonlinear vibration by using the Galerkin’s and multiple scale methods for the FGM thermo‐ electro‐elastic nanostructures. In 2021, Guo et al. [5] presented the non-dimensional frequency by using the element-free IMLS-Ritz method for the graphene nanoplatelets reinforced composite (GPLRC) plates with matrix cracks. There are some non-dimensional and dimensional frequency of vibration presentations in the composite structure of plates. In 2020, Babaei et al. [6] used the 3D elasticity theory of structures to investigate the non-dimensional natural frequency of saturated porous annular sector in thick FGM plates. Some numerical solution of finite element method (FEM) are obtained. In 2019, Alaimo et al. [7] used the analytical procedure of Navier to investigate the frequency of composite plates embedding viscoelastic layers. Some numerical results of FEM are obtained. In 2019, Safaei et al. [8] used the mesh-free method and HSDT to investigate the frequency of FGM carbon nanotube (CNT) plate under thermo-mechanical loads. Some numerical results of dynamic behaviors are obtained. In 2019, Sun and Wei [9] used the singular boundary method (SBM) to investigate the frequency of thin elastic plates. Some numerical results of harmonic behaviors are obtained. In 2019, Vinyas et al. [10] used the HSDT of displacements to investigate the frequency of three-phase smart magneto-electro-elastic (TPS-MEE) composite plates. Some numerical results of FEM are obtained. In 2017, Rezaei et al. [11] investigated the non-dimensional frequency of * Corresponding author. E-mail address: cchong@hust.edu.tw Tel.: +886-919037599; Fax: +886-4-24961187 Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 22 FGM plates with porosities by using linear model of displacements with analytical method. Some numerical results of non-dimensional frequencies are obtained. In 2017, Yao et al. [12] used von Karman large deformation theory and Hooke law to investigate the frequency responses of the bistable piezoelectric plate. Some numerical results of frequencies are obtained. Jha et al. investigated the numerical solutions of thick FGM with the higher order shear/shear-normal deformations theories (HOSTs/HOSNTs) in 2013, the numerical frequency solutions for free vibration of simply supported FGM plates are obtained [13]. In 2017, Duc et al. [14] used the FSDT displacement method to present the analytical solutions for static response and free vibration of carbon nanotube reinforced composite (CNTRC) FGM plates. Several numerical investigations have been studied in the thick and thin FGMs. Hong [15] used the nonlinear model TSDT to obtain the vibration solutions for the thermal stress and deflection of the thick FGM plates in 2019. It would be interesting to study the non-dimensional and dimensional natural frequency containing the c1 displacement term in the TSDT displacement model and calculated values of shear correction for the thick FGM plates under free vibration. Length to thickness ratio, temperature of environment and power law index of FGM three parametric effects on the non-dimensional and dimensional natural frequency of thick FGM plates are also studied. 2. Procedure of Formulations The typically two-material plates respectively in constituent material 1 and constituent material 2 of FGMs are shown in Fig. 1 with thickness h1 and h2 are considered for the derivation of formulation under the thermal influence of environment. Usually, power-law type of materials function are selected to study the numerical solutions for the two-material. These standard forms of properties e.g. Young’s modulus of FGM Efgm can be denoted by Efgm = (E2-E1)[(z+h * /2)/h * ] Rn , in which Rn is power index, E1 and E2 are Young’s modulus of constituent material 1 and 2, respectively, for the other properties can be simplified into the mean expressions e.g. Poisson’s ratios of FGM is νfgm = (ν1+ν2)/2, in which ν1 and ν2 are Poisson’s ratios of constituent material 1 and 2, respectively [16]. Also the constituent material properties Pi are denoted and found in the variation of temperature T of environment i.e. Pi = P0(P-1T -1 +1+P1T+P2T 2 +P3T 3 ), in which P0, P-1, P1, P2 and P3 are temperature coefficients. Displacements in thick plates with variables u, v and w of FGM are in dependent of time. These displacements can be denoted and assumed in the additional nonlinear coefficient c1 term for c1 = 4/[3(h * ) 2 ] presented in TSDT mode [17], where h * is the summation of two-material thickness. Fig. 1 Two-material FGM plate under thermal environment The nonlinear Reddy’s TSDT displacement field, strain-displacement relations and dynamic equilibrium differential equations of motion governing the FGM plates under consideration of formulations can be referred to the published paper in 2019 by Hong [15]. The inter-laminar stresses within thick FGM plate under the influence of temperature difference △T for the k th constituent layer can be derived and given [18-19]. The TSDT term of displacement modes in thick FGM plate contained in the dynamic equilibrium equations can be derived and given [20] with k th constituent ply density ρ(k) integration parameters Ii, Ji and K2, in the subscript i = 0, 1, 2,…, 6. Von Karman type of strain-displacement relation equations can be derived and used with displacements u 0 and v 0 are in the axes direction x and y respectively, with transverse displacement w is in the axis direction z within the middle plate of thick FGM, also with the shear rotations Ψx and Ψy are in the x and y direction of axes respectively. By substituting stress and Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 23 strain-displacement equations into dynamic equations of motion, thus the dynamic differential equations of equilibrium contained TSDT terms of thick FGM plates can be obtained in matrix forms in terms of partial derivatives of displacements, shear rotations, mechanical loads (p1 and p2 are external loads in axes direction x and y respectively, q is pressure load), thermal loads and inertia terms. These calculated correction values of shear coefficient kαhare typically only in functions of parameters T, Rn and h * [15]. Stiffness parameters are integrated and given as follows. * * / 2 2 3 4 6 / 2 ( , , , , , ) (1, , , , , ) s s s s s s s s s s s s s s h i j i j i j i j i j i j h i j A B D E F H Q z z z z z dz − = ∫ (1) * * * * * * * * * * * * * * * * / 2 2 3 4 / 2 5 ( , , , , , ) (1, , , , , ) h i j i j i j i j i j i j h i j A B D E F H k Q z z z z z dz − = ∫ α (2) where ������ and ���∗�∗ are transformed reduced stiffness of FGM can be simplified in 2007 by Shen [21], in which i s , j s = 1, 2, 6 and i * , j * = 4, 5. 3. Numerical Results and Discussion The two-material FGM thick plates with same direction of constituents are selected to investigate the non-dimensional and dimensional frequency of numerical results under free vibration without thermal forces, distributed forces and pressure forces (△T = 0, p1 = p2 = 0, q = 0). The dimensional frequency ωmn (the subscripts m and n are mode shape numbers in the x and y direction of axes) can be simplified (I1 = I3 = J1 = 0, Bij = Eij = 0, D16 = D26 = A16 = A26 = 0 and A45 = D45 = F45 = 0) and derived for the following displacements and shear rotations in simply supported boundary equations, 0 ( / ) ( / ) ( ) mn mn u a cos m x a sin n y b sin t= π π ω (3) 0 ( / ) ( / ) ( ) mn mn b sin m x a cos n y b sin tv = π π ω (4) ( / ) ( / ) ( ) mn mn c sin m x a sin n y b sinw t= π π ω (5) ( / ) ( / ) ( ) x mn mn d s m x a sin n y b sio n tc=ψ π π ω (6) ( / ) ( / ) ( ) y mn mn e sin m x a cos n y b sin t=ψ π π ω (7) in which amn, bmn, cmn, dmn, and emn are amplitudes, t is the time, a is the length and b is the width of thick FGM plates. The major difference in the present free vibration work of simply homogeneous equation with respect to fully homogeneous equation is the calculated simplification in the matrix element terms. After Eqs. (3)-(7) have been substituted into equations of equilibrium with no external loads (f1 = f2 =...= f5 = 0), and by conveniently assuming that matrix element terms FH13 = FH14 = FH15 = FH23 = FH24 = FH25 = 0 and I3 = J1 = I6 = J4 = 0 in the fully homogeneous equation, thus the simply homogeneous equation for the dimensional frequency ωmn contained the c1 terms of TSDT can be obtained as follows [15]. 11 12 12 22 33 34 35 2 34 44 45 0 2 35 45 55 0 0 0 0 00 0 0 00 0 0 0 0 0 0 0 0 mn mnmn mnmn mn mn mn mn mn FH FH aFH FH bFH FH FH cK FH FH FH I d eK FH FH FH I − − − = − −                                                  λ λ λ λ λ (8a) Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 24 where 0 2 mn mn I=λ ω (8b) 2 2 11 11 66 ( / ) ( / )FH A m a A n b= +π π (8c) 1 612 62 )( ( / )( / )FH A A m a n b= + π π (8d) 22 66 2 2 2 2 ( / ) ( / )FH A m a A n b= +π π (8e) 2 2 2 4 2 2 2 2 33 55 44 1 11 1 1 66 2 4 2 2 22 1 55 1 55 1 44 1 44 ( ) ( ) ( / ) (2 4 )( / ) ( / ) ( / ) 3 (2 3 / / )( / ) 3 (2 3 )( / ) FH A m a A n b c H m a c c H m a n b c H n b c D c F m a c D c F n b = + + + + + − − − − π π π π π π π π (8f) In the derivation of mathematics, the algebraic determinant of Eq. (8a) could be presented in the following equation in terms of λmn, thus the ωmn can be calculated. 5 4 3 2 (1) (2) (3) (4) (5) (6) 0 mn mn mn mn mn A A A A A A+ + + + + =λ λ λ λ λ (9a) where (1)A sd= − (9b) 11 22 (2) ( )A FH FH sd sc+= + (9c) 22 1211 1212 11 (3) [( ) ( ) ]A FH FH FH FH sd FH FH sc sb− −= ++ + (9d) 22 111 12 1 212 2 (4) ( ) ( )A FH FH FH FH sc FH FH sb sa−= + + + (9e) 22 1211 22 21 11 (5) [( ) ( ) ]A FH FH FH FH sb FH FH sa+− += − (9f) 22 11 12 21 (6) ( )A FH FH FH FH sa−= (9g) in which 2 2 0 ( / )sd K I= (9h) 33 44 2 0 /sdsc FH FH K I+= (9i) 33 55 44 55 33 44 35 35 34 34 452 0 45 ( /)sb FH FH FH FH FH FH FH FH FH FH K I FH FH= −+ + − − (9j) 33 44 55 44 34 35 35 34 45 35 35 44 34 34 55 45 45 33 sa FH FH FH FH FH FH FH FH FH FH FH FH FH FH FH FH FH FH + + − −= − (9k) A numerical method used to calculate and solve the polynomial Eq. (9a) in terms of fifth-order of � , thus the natural frequency ωmn can be determined. A FORTRAN program coded in the algorithm of Newton's method by Conte and de Boor in 1980 [22] is used to find the solution of Eq. (9a). The accuracy of tolerance is 1e-06 and used in the program to meet the condition of convergence then find the root of Eq. (9a). The composited FGM two-material SUS304/Si3N4 is implemented for the computation results under temperature of environment T. The material 1 SUS304 (stainless steel) is located at lower Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 25 position, the material 2 Si3N4 (silicon nitride) is located at upper position of the FGM thickness, respectively. In the preliminary and fundamental calculations, it did not consider the c1 effect of TSDT onto the derivative of shear coefficient kα. For example a/b = 1, h1 = h2, h * = 0.12mm-12mm, the kα values under T = 300K, 600K and 1000K are shown in Table 1(a), Table 1(b) and Table 1(c), respectively. When the values of kα are calculated under the condition of cases only for the corresponding T, Rn and h * , then they are used in getting results for Figs. 2-5. Table 1(a) Varied kα vs. Rn and h * under T = 300K c1(1/mm 2 ) h * (mm) kα Rn: 0.1 0.2 0.5 1 2 5 10 92.592598 0.12 0.098323 0.127625 0.266187 0.638005 1.007910 0.932667 0.843506 0.925925 1.2 0.079956 0.084882 0.102678 0.138574 0.217517 0.396486 0.492255 0.231481 2.4 0.075063 0.074789 0.074832 0.074112 0.064424 0.021065 0.001152 0.037037 6 0.069011 0.063127 0.048652 0.030937 0.010870 0.000220 0.000000 0.009259 12 0.064733 0.055449 0.034895 0.015697 0.002738 0.000006 0.000000 Table 1(b) Varied kα vs. Rn and h * under T = 600K c1(1/mm 2 ) h * (mm) kα Rn: 0.1 0.2 0.5 1 2 5 10 92.592598 0.12 0.085699 0.111597 0.236756 0.596251 0.995481 0.923808 0.834206 0.925925 1.2 0.069478 0.073722 0.089181 0.120708 0.191371 0.359433 0.456324 0.231481 2.4 0.065173 0.064854 0.064737 0.063996 0.055506 0.018026 0.000982 0.037037 6 0.059859 0.054641 0.041932 0.026560 0.009297 0.000188 0.000000 0.009259 12 0.056109 0.047937 0.030016 0.012979 0.002340 0.000005 0.000000 Table 1(c) Varied kα vs. Rn and h * under T = 1000K c1(1/mm 2 ) h * (mm) kα Rn: 0.1 0.2 0.5 1 2 5 10 92.592598 0.12 0.071316 0.091135 0.192482 0.562116 1.236624 1.094820 0.917769 0.925925 1.2 0.057519 0.059382 0.067907 0.088108 0.137812 0.270674 0.355043 0.231481 2.4 0.053882 0.052072 0.048781 0.045488 0.037321 0.011122 0.000565 0.037037 6 0.049404 0.043712 0.031294 0.018557 0.006113 0.000115 0.000000 0.009259 12 0.046254 0.038257 0.022289 0.009340 0.001533 0.000003 0.000000 The non-dimensional frequency parameters f * , ω*, Ω in terms of fundamental dimensional frequency ω11 with mode shape numbers in the subscripts m = n = 1 are used for SUS304/Si3N4 thick plate under free vibration and given as follows, 2 * * 11 2 /f h E= ρω (10) * 1 2 2 1 / ) /( s s b I D= ωω π (11) * 2 1 111 2 1 / ) (1 /( )a h E−Ω = ρ νω (12) * * / 2 1 / 2 h h s I dz − = ∫ ρ (13) * * / 2 2 1 / 2 h s h D z dzQ − = ∫ (14) 2 1 1 1 / (1 )Q E= −ν (15) in which ω11 is the fundamental dimensional frequency, ρ2 and ρ1 are the density of FGM material 2 and material 1, respectively. In the Tables 2-4, the present results of non-dimensional f * , ω* and Ω are compared with available published paper. In Table 2, the present results of f * vs. h * are shown for a/h * =10, 300K, kα and c1 values, we found f * = 0.080832 at h * = 900mm, Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 26 Rn = 0.5 is in close to f * = 0.0839 for Al/ZrO2 by Jha et al. in 2013 [13] with HOSNT12 in which 12 degrees of freedom considered for the model. In Table 3, the present results of ω* vs. h * are shown for a/h * =10, 300K, kα and c1 values, we found ω * = 4.094445 at h * = 600mm, Rn = 2 is in close to ω * = 4.1165 with h * = 200mm forced vibration by Kim in 2005 [23]. Table 2 Comparison of non-dimensional f * c1(1/mm 2 ) h * (mm) f * Present results, a/ h * = 10, T = 300K, varied kα, for SUS304/Si3N4 Jha et al. 2013, for Al/ZrO2, Rn = 0.5 Rn = 0.5 Rn = 1 Rn = 2 0.925925 1.2 0.003713 0.002326 0.002412 - 0.333333 2 0.008001 0.004982 0.005163 - 0.000033 200 0.017724 0.017726 0.017708 - 0.000014 300 0.026743 0.026744 0.026727 - 0.000003 600 0.053789 0.053789 0.053773 - 0.000001 900 0.080832 0.080832 0.080816 0.0839 Table 3 Comparison of non-dimensional ω* c1(1/mm 2 ) h * (mm) ω* Present results, a/ h * = 10, T = 300K, varied kα Kim 2005, Forced vibration, h * = 200mm, Rn = 2 Duc et al. 2017 CNTRC, FSDT Rn: 0.5 1 2 0.925925 1.2 0.282727 0.177156 0.183727 - - 0.333333 2 0.609222 0.379389 0.393123 - - 0.000033 200 1.349598 1.349716 1.348404 - - 0.000014 300 2.036305 2.036380 2.035098 - - 0.000003 600 4.095676 4.095697 4.094445 4.1165 3.99244 (UD) 0.000001 900 6.154796 6.154793 6.153549 - - Also, compare the present results with the analytical FSDT results of uniform distribution (UD) in CNTRC FGM plates by Duc et al. in 2017 [14]. In Table 4, the present results of Ω vs. h * are shown for a/h * = 10, 700K, kα and c1 values, we found Ω = 5.1314993 at h * = 250mm, Rn = 1 is in close to Ω = 5.359 with forced vibration under temperature rise (△T = 400K) by Ungbhakorn & Wattanasakulpong in 2013 [24]. Table 4 Non-dimensional Ω comparison for SUS304/Si3N4 c1(1/mm 2 ) h * (mm) Ω Present results, a/h * = 10, T = 700K, varied kα Ungbhakorn & Wattanasakulpong 2013, Rn = 1 Rn: 0.5 1 2 0.925925 1.2 0.548323 0.568358 0.983723 - 0.333333 2 1.966345 1.217748 1.266043 - 0.000033 200 4.089171 4.089357 4.084728 - 0.000021 250 5.131368 5.131493 5.126937 5.359 0.000014 300 6.173302 6.173382 6.168875 - The dimensional ωmn (unit 1/s) under no temperature difference (△T = 0) of two directional mode in the subscripts m, n are investigated. In Fig. 2, the present results of ω1n vs. Rn are shown for a/h * = 5 and 10, respectively for T = 300K, kα and c1 values. For the cases of Rn = 0.5, 1 and 10, the present results of ω1n are oscillating and converging to nearly 0.005/s with n (in the subscript n = 1-9). The greatest value of ω11 = 0.054676/s is found, then decreasing to value ω19 = 0.005694/s for a/h * = 5, Rn = 10. The greatest value of ω11 = 0.067990/s is found, then decreasing to value ω19 = 0.005715/s for a/h * = 10, Rn = 10. When n = 9, the frequencies are almost meeting at a point, e.g. ω19 = 0.005694/s for a/h * = 5, Rn = 0.5, 1 and 10, also ω19 = 0.005715/s for a/h * = 10, Rn = 0.5 and 1, but not for Rn = 10. The frequencies would be meeting closely at a point for all the value of Rn when the values of n are greater more than 9. Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 27 (a) a/h * =5 (b) a/h * =10 Fig. 2 ω1n (unit 1/s) vs. Rn In Fig. 3, the present results of ω1n vs. T are shown for a/h * = 5 and 10, respectively for Rn = 0.5, kα and c1 values. Generally the values of ω1n are oscillating and converging to nearly 0.005/s with n (in the subscript n = 1-9) for T = 300K, 600K and 1000K. The greatest value of ω12 = 0.059699/s is found, then decreasing to value ω19 = 0.001000/s for a/h * = 5, T = 600K. The greatest value of ω11 = 0.060500/s is found, then decreasing to value ω19 = 0.010078/s for a/h * = 10, T = 600K. When n = 9, the frequencies are almost meeting at a point, e.g. ω19 = 0.001000/s for a/h * = 5, T = 600K and 1000K, but not for T = 300K, also ω19 = 0.010078/s for a/h * = 10, T = 300K, 600K and 1000K. The frequencies would be meeting closely at a point for all the value of T when the values of n are greater more than 9. (a) a/h * = 5 (b) a/h * = 10 Fig. 3 ω1n (unit 1/s) vs. T The influence of gradation of material properties e.g. Rn values and thermal effect e.g. T values used on the vibration can be described in Figs. 2-3 as follows. The present results of ω1n are oscillating and converging with the Rn values. Also the present results of ω1n are oscillating and converging with the T values. The trend of the curves will not change for the non-dimensional and dimensional frequencies. Fig. 4 ω9n (unit 1/s) vs. n Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 28 In Fig. 4, the present results of ω9n vs. n are shown for a/h * = 5 and 10, respectively for T = 300K, Rn = 0.5, kα and c1 values. In general, the values ω9n are oscillating with n (in the subscript n = 1-9) for a/h * = 10. The values ω9n are linearly decreasing with n for a/h * = 5. The greatest value of ω94 = 0.011816/s is found, then decreasing to value ω97 = 0.001000/s for a/h * = 10. The greatest value of ω91 = 0.005274/s is found, then decreasing to value ω99 = 0.003356/s for a/h * = 5. For TSDT coefficient c1 = 0, the results can be added as shown in Fig. 5 and recently can't be compared with pretty rare published data. In Fig. 5(a), the present results of ω1n vs. Rn are shown for a/h * = 5, T=300K, varied kα and c1 = 0. The values of ω1n are decreasing smoothly with n = 1-6, but have a jump at n = 7, then also decreasing smoothly with n = 8-9. The greatest value of ω1n is 0.019021/s for Rn = 10 at n = 1, it is almost three times smaller than that 0.054676/s in the nonzero of c1 case. In Fig. 5(b), the present results of ω1n vs. T are shown for a/h * = 5, Rn = 0.5, varied kα and c1 = 0. The values of ω1n are decreasing smoothly with n = 1-6 for T = 300K and 600K, but have a jump at n = 7, then also decreasing smoothly with n = 8-9. The values of ω1n are decreasing smoothly with n = 1-5 for T = 1000K, but have a jump at n = 6, then also decreasing smoothly with n = 7-9. The greatest value of ω1n is 0.0188063/s for T = 1000K at n = 1, it is a little smaller than that 0.020803/s in the nonzero of c1 case. In Fig. 5(c), the present results of ω9n vs. n are shown for a/h * = 5, T = 300K, Rn = 0.5, varied kα and c1 = 0. The values of ω9n are decreasing linearly with n = 1-9. The greatest value of ω9n is 0.007045/s at n = 1, it is a little greater than that 0.005274/s in the nonzero of c1 case. The varied kα and c1 value of TSDT have great effect on the values of natural frequencies. (a) ω1n (unit 1/s) vs. Rn for T = 300K (b) ω1n (unit 1/s) vs. T for Rn = 0.5 (c) ω9n (unit 1/s) vs. n for T=300K, Rn = 0.5 Fig. 5 c1=0 for a/h * = 5 and varied kα Similar way for keeping n = 1 and varying m from 1 to 9, the ωm1 can also be plotted and results can be included and discussed for varied kα and c1 values as shown in Figs. 6-7. In Fig. 6, the present results of ωm1 vs. Rn are shown for a/h * = 5 and 10, respectively for T = 300K, kα and c1 values. The cases of Rn = 0.5, 1 and 10 in the present results of ωm1 are small oscillating and converging to nearly 0.005/s with m (in the subscript m = 1-9). The greatest value of ω11 = 0.054676/s is found, then decreasing to value ω91 = 0.005969/s for a/h * = 5, Rn = 10. The greatest value of ω11 = 0.067990/s is found, then decreasing to value ω91 = 0.004257/s for a/h * = 10, Rn = 10. When m = 9, the frequencies are almost meeting at a point, e.g. ω91 = 0.005969/s Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 29 for a/h * = 5, Rn = 0.5, 1 and 10, also ω91 = 0.006950/s for a/h * = 10, Rn = 0.5 and 1, but not for Rn = 10. The frequencies would be meeting closely at a point for all the value of Rn when the values of m are greater more than 9. In Fig. 7, the present results of ωm1 vs. T are shown for a/h * = 5 and 10, respectively for Rn = 0.5, kα and c1 values. The cases of T = 300K, 600K and 1000K in the present results of ωm1 are small oscillating and converging to nearly 0.005/s with m (in the subscript m = 1-9). The greatest value of ω11 = 0.031320/s is found, then decreasing to value ω91 = 0.005513/s for a/h * = 5, T = 600K. The greatest value of ω11 = 0.098299/s is found, then decreasing to value ω91 = 0.007040/s for a/h * = 10, T = 1000K. When m = 9, the frequencies are almost meeting at a point, e.g. ω91 = 0.005513/s for a/h * = 5, T = 300K and 600K, but not for T = 1000K, also ω91 = 0.007040/s for a/h * = 10, T = 300K, 600K and 1000K. The frequencies would be meeting closely at a point for all the value of T when the values of m are greater more than 9. (a) a/h * = 5 (b) a/h * = 10 Fig. 6 ωm1 (unit 1/s) vs. Rn (a) a/h * = 5 (b) a/h * = 10 Fig. 7 ωm1 (unit 1/s) vs. T The results of non-dimensional and dimensional natural frequency are obtained in free vibration of two-material thick FGM plates with power law index by considering the TSDT model, varied shear coefficient and temperature of environment. Some of the present results of non-dimensional f * , ω* and Ω are compared with available published paper and found in close data. It is fundamental to study the effects of varied shear coefficient values and additional c1 displacement term in nonlinear coefficient of TSDT on the frequency calculations under free vibration. The varied kα and c1 value of TSDT have great effect on the values of natural frequencies. For the cases of Rn = 0.5, 1 and 10, the present results ω1n are oscillating and converging to nearly 0.005/s with n (in the subscript n = 1-9). For the cases of T = 300K, 600K and 1000K, the present results ω1n are oscillating and converging to nearly 0.005/s with n (in the subscript n = 1-9). Similar way for keeping n = 1 and varying m from 1 to 9, the cases of Rn = 0.5, 1 and 10 in the results of ωm1 are small oscillating and converging to nearly 0.005/s with m (in the subscript m = 1-9), the cases of T = 300K, 600K and 1000K in the results of ωm1 are small oscillating and converging to nearly 0.005/s with m (in the subscript m = 1-9). Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 30 4. Conclusions The natural frequency of free vibration in two-material thick FGM plates are great effected by power law index, the TSDT model, varied shear coefficient and temperature of environment. It is fundamental and new to study the additional effect of nonlinear displacement TSDT c1 term on the frequency for two-material thick FGM plates. The present natural frequency results are oscillating and converging with respect to mode shapes by simultaneously including x and y directions. These natural frequency data of free vibration can provide a base reference and used further in the thermal vibration study in the future. Nomenclature h1 Thickness of constituent material 1 h2 Thickness of constituent material 2 h * Summation of two-material thickness a Length of thick FGM plates b Width of thick FGM plates x x direction of Cartesian axes y y direction of Cartesian axes z z direction of Cartesian axes T Temperature of environment Rn Power index in power-law type of FGM materials function Efgm Young’s modulus of FGM E1 Young’s modulus of constituent material 1 E2 Young’s modulus of constituent material 2 νfgm Poisson’s ratios of FGM ν1 Poisson’s ratios of constituent material 1 ν2 Poisson’s ratios of constituent material 2 ρ1 Density of constituent material 1 ρ2 Density of constituent material 2 Pi Constituent material properties P0, P-1, P1, P2 and P3 Temperature coefficients of constituent materials u, v and w Displacements in x, y and z directions of Cartesian axes u 0 and v 0 Displacements in the axis directions x and y within the middle plate of thick FGM c1 Nonlinear Reddy’s TSDT displacement field c1 = 4/[3(h * ) 2 ] △T Temperature difference k th k th constituent layer ρ(k) kth constituent ply density Ii, Ji and K2 Density ρ (k) integration parameters, in the subscript I = 0, 1, 2,…, 6 Ψx and Ψy Shear rotations in the x and y direction of axes p1 and p2 External loads in axes direction x and y q Pressure load kα Shear coefficient ������ and ���∗�∗ Transformed reduced stiffness of FGM with i s , j s = 1, 2, 6 and i * , j * = 4, 5 ωmn Dimensional frequency with the subscripts m and n are mode shape numbers in the x and y direction of axes ω11 Fundamental dimensional frequency with m = n = 1 f * , ω*, Ω Non-dimensional frequency parameters t Time f1, f2,..., f5 External loads λmn Frequencies parameter in λmn = I0ωmn 2 Conflicts of Interest The authors declare no conflict of interest. Proceedings of Engineering and Technology Innovation, vol. 17, 2021, pp. 21-31 31 References [1] K. K. Żur, M. Arefi, J. Kim, and J. N. Reddy, “Free Vibration and Buckling Analyses of Magneto-Electro-Elastic FGM Nanoplates Based on Nonlocal Modified Higher-Order Sinusoidal Shear Deformation Theory,” Composites Part B, vol. 182, 107601, February 2020. [2] H. M. Ouakad, A. Valipour, K. K. Żur, H. M. Sedighi, and J. N. Reddy, “On the Nonlinear Vibration and Static Deflection Problems of Actuated Hybrid Nanotubes Based on the Stress-Driven Nonlocal Integral Elasticity,” Mechanics of Materials, vol. 148, 103532, September 2020. [3] A. Ghobadi, H. Golestanian, Y. T. Beni, and K. K. Żur, “On the Size-Dependent Nonlinear Thermo-Electro-Mechanical Free Vibration Analysis of Functionally Graded Flexoelectric Nano-Plate,” Communications in Nonlinear Science and Numerical Simulation, in press. [4] A. Ghobadi, Y. T. Beni, and K. K. Żur, “Porosity Distribution Effect on Stress, Electric Field and Nonlinear Vibration of Functionally Graded Nanostructures with Direct and Inverse Flexoelectric Phenomenon,” Composite Structures, in press. [5] H. Guo, T. Yang, K. K. Żur, and J. N. Reddy, “On the Flutter of Matrix Cracked Laminated Composite Plates Reinforced with Graphene Nanoplatelets,” Thin-Walled Structures, vol. 158, 107161, January 2021. [6] M. Babaei, M. H. Hajmohammad, and K. Asemi, “Natural Frequency and Dynamic Analyses of Functionally Graded Saturated Porous Annular Sector Plate and Cylindrical Panel Based on 3D Elasticity,” Aerospace Science and Technology, vol. 96, 105524, January 2020. [7] A. Alaimo, C. Orlando, and S. Valvano, “Analytical Frequency Response Solution for Composite Plates Embedding Viscoelastic Layers,” Aerospace Science and Technology, vol. 92, pp. 429-445, September 2019. [8] B. Safaei, R. Moradi-Dastjerdi, Z. Qin, and F. Chu, “Frequency-Dependent Forced Vibration Analysis of Nanocomposite Sandwich Plate under Thermo-Mechanical Loads,” Composites Part B, vol. 161, pp. 44-54, March 2019. [9] L. Sun and X. Wei, “A Frequency Domain Formulation of the Singular Boundary Method for Dynamic Analysis of Thin Elastic Plate,” Engineering Analysis with Boundary Elements, vol. 98, pp. 77-87, January 2019. [10] M. Vinyas, K. K. Sunny, D. Harursampath, T. Nguyen-Thoi, and M. A. R. Loja, “Influence of Interphase on the Multi-Physics Coupled Frequency of Three-Phase Smart Magneto-Electro-Elastic Composite Plates,” Composite Structures, vol. 226, 111254, October 2019. [11] A. S. Rezaei, A. R. Saidi, M. Abrishamdari, and M. H. Pour Mohammadi, “Natural Frequencies of Functionally Graded Plates with Porosities via a Simple Four Variable Plate Theory: An Analytical Approach,” Thin-Walled Structures, vol. 120, pp. 366-377, November 2017. [12] M. Yao, W. Hu, and W. Zhang, “Nonlinear Frequency Responses of the Bistable Piezoelectric Plate,” Procedia IUTAM, vol. 22, pp. 208-215, 2017. [13] D. K. Jha, T. Kant, and R. K. Singh, “Free Vibration Response of Functionally Graded Thick Plates with Shear and Normal Deformations Effects,” Composite Structures, vol. 96, pp. 799-823, February 2013. [14] N. D. Duc, J. Lee, T. Nguyen-Thoi, and P. T. Thang, “Static Response and Free Vibration of Functionally Graded Carbon Nanotube-Reinforced Composite Rectangular Plates Resting on Winkler–Pasternak Elastic Foundations,” Aerospace Science and Technology, vol. 68, pp. 391-402, September 2017. [15] C. C. Hong, “GDQ Computation for Thermal Vibration of Thick FGM Plates by Using Fully Homogeneous Equation and TSDT,” Thin-Walled Structures, vol. 135, pp. 78-88, February 2019. [16] S. H. Chi and Y. L. Chung, “Mechanical Behavior of Functionally Graded Material Plates Under Transverse Load, Part I: Analysis,” International Journal of Solids and Structures, vol. 43, no. 13, pp. 3657-3674, June 2006. [17] S. J. Lee, J. N. Reddy, and F. Rostam-Abadi, “Transient Analysis of Laminated Composite Plates with Embedded Smart-Material Layers,” Finite Elements in Analysis and Design, vol. 40, no. 5-6, pp. 463-483, March 2004. [18] S. J. Lee and J. N. Reddy, “Non-Linear Response of Laminated Composite Plates Under Thermomechanical Loading,” International Journal of Non-Linear Mechanics, vol. 40, no. 7, pp. 971-985, September 2005. [19] J. M. Whitney, “Structural Analysis of Laminated Anisotropic Plates,” U.S.A: CRC Press, 1987. [20] J. N. Reddy, “Energy Principles and Variational Methods in Applied Mechanics,” New York: Wiley, 2017. [21] H. S. Shen, “Nonlinear Thermal Bending Response of FGM Plates due to Heat Condition,” Composites Part B: Engineering, vol. 38, pp. 201-215, 2007. [22] S. D. Conte and C. de Boor, “Elementary Numerical Analysis, an Algorithmic Approach,” 3rd ed. U.S.A.: McGraw-Hill Book Company, 1980. [23] Y. W. Kim, “Temperature Dependent Vibration Analysis of Functionally Graded Rectangular Plates,” Journal of Sound and Vibration, vol. 284, no. 3-5, pp. 531-549, June 2005. [24] V. Ungbhakorn and N. Wattanasakulpong, “Thermo-Elastic Vibration Analysis of Third-Order Shear Deformable Functionally Graded Plates with Distributed Patch Mass Under Thermal Environment,” Applied Acoustics, vol. 74, no. 9, pp. 1045-1059, September 2013. Copyright© by the authors. Licensee TAETI, Taiwan. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY-NC) license (https://creativecommons.org/licenses/by-nc/4.0/).