Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 56, 4, pp. 1163-1178, Warsaw 2018 DOI: 10.15632/jtam-pl.56.4.1163 BENDING, BUCKLING AND FREE VIBRATION OF A BEAM WITH UNSYMMETRICALLY VARYING MECHANICAL PROPERTIES Krzysztof Magnucki, Jerzy Lewiński Institute of Rail Vehicles TABOR, Poznań, Poland Ewa Magnucka-Blandzi Poznan University of Technology, Institute of Mathematics, Poznań, Poland e-mail: ewa.magnucka-blandzi@put.poznan.pl Piotr Kędzia Poznan University of Technology, Institute of Applied Mechanics, Poznań, Poland The subject of the paper is a beamwith unsymmetrically varying mechanical properties in the depth direction.The nonlinear hypothesis of plane cross section deformation is assumed. BasedonHamilton’s principle, twodifferential equations ofmotion are obtained.The system of equations is analytically solved with a view to analyse the bending, buckling and free vibration problems of the beam. Moreover, the FEM model of the beam is developed and deflections, critical axial forces and natural frequencies of the beam are calculated. The results of these twomethods are compared. Keywords: FGMbeam, mathematical modelling, numerical FEM calculations 1. Introduction Elementswithvaryingmechanical properties areapplied inmodernconstructions.Kubiak(2005) presented dynamic buckling problems of thin-walled composite plates with varying width-wise material properties. Zhang et al. (2006) presented free vibration analysis of rectangular com- posite laminated plates. Zenkour (2006) analysed bending problems of rectangular functionally graded plates under a transverse uniform load. Birman and Byrd (2007) presented a review of the papers published since 2000 related to the modelling and analysis of functionally graded materials and structures. Kapuria et al. (2008) described the theoretical model of bending and free vibration of layered functionally graded beams and its experimental validation. Debowski et al. (2010) studied the dynamic stability problem of a metal foam rectangular plate under com- pression in themiddle plane.Magnucka-Blandzi (2011) presented bending anddynamic stability results of studies of the sandwich beam with a metal foam core. Kubiak (2011) described an estimation problem of dynamic buckling for composite columns with open cross-sections. Thai and Vo (2012, 2013) presented bending, buckling, and vibration of functionally graded beams and plates with the use of nonlinear shear deformation theories. Mahi et al. (2015) presented bending and free vibration analysis of isotropic, functionally graded sandwich and laminated composite plates with the use of a new hyperbolic deformation theory. Kolakowski and Mania (2015) presented the dynamic response of thin functionally graded plates with a static unsym- metrical stable postbuckling path. Chen et al. (2015, 2016a,b) analysed static bending, elastic buckling and free vibrations problems of shear deformable functionally graded porous beams and sandwich beams with a functionally graded porous core. Jun et al. (2016) studied the free vibration problem of axially loaded laminated composite beams using a unified higher-order shear deformation theory and a dynamic stiffness method. Mojahedin et al. (2016) presented 1164 K.Magnucki et al. the buckling problem of functionally graded circular plates with symmetrically and unsymme- trically varyingmechanical properties based on a higher order shear deformation theory. Li and Hu (2016) analysed nonlinear bending and free vibration problems of nonlocal strain gradient beams made of a functionally graded material. Feyzi and Khorshidvand (2017) presented the axisymmetric post-buckling behaviour problem of saturated porous circular plates. Song et al. (2017) described vibration problems of functionally graded polymer composite plates reinforced with graphene nanoplatelets. Smyczynski andMagnucka-Blandzi (2018) presented a comparison of the study results of three-point bending of a sandwich beam with two binding layers with the use of two nonlinear hypotheses. Sayyad and Ghugal (2017) presented an extensive review of the papers devoted to bending, buckling and free vibration problems with special attention paid to the shear effects. The subject of the study is a beam with unsymmetrically varying mechanical properties. A nonlinear hypothesis of deformation of the plane cross section of the beam is developed. Particular attention is paid to location of the neutral axis with consideration of the shear effect. Variability of the elasticmodulus –Young’smodulus in the depthdirection of the beam is shown in Fig. 1. Fig. 1. Scheme of the elastic modulus variability in the depth direction of the beam The values of elasticity moduli andmass density of the beam vary as follows E(y) = 1 2 E11+e2− (1−e2)sin(πη)] G(y) = 1 2 G1[1+g2− (1−g2)sin(πη)] ρ(y)= 1 2 ρ1[1+ ρ̃2− (1− ρ̃2)sin(πη)] (1.1) where: e2 = E2/E1, g2 = G2/G1 = (1+ ν1)/(1+ ν2)e2, ρ̃2 = ρ2/ρ1 are dimensionless relative parameters, E1, E2 – Young’s moduli, ν1, ν2 – Poisson’s ratios, ρ1, ρ2 –mass densities, η = y/h – dimensionless coordinate (−0.5¬ η ¬ 0.5), h – depth of the beam. The relationship between the relative density andYoung’smoduli ratio ρ̃2 = √ e2 is assumed based on the papers by Chen et al. (2015, 2016b). 2. Analytical model of the beam The nonlinear hypothesis is assumed for the purpose of modelling of the beam. A plane cross section before bending is no longer plane after bending of the beam (Fig. 2). This hypothesis is a generalization of the shear deformation theory for functionally graded structures. Two coordinate systems are adopted – x,y and x1,y1 (Fig. 2). The x1 axis is the neutral axis, therefore, thedisplacement v(x1, t) is equivalentwith v(x,t).The coordinate y1 = h(η+η0), Bending, buckling and free vibration of a beam... 1165 Fig. 2. Deformation of the plane cross section of the beam – the nonlinear hypothesis where η0 = y0/h, therefore, based on the above hypothesis, the displacement is in the following form u(x,y,t) =−h { (η+η0) ∂v ∂x − [sin(πη)+sin(πη0)]ψ(x,t) } (2.1) where: v(x,t) – deflection, ψ(x,t) – dimensionless function of the shear effect. The shear effect displacements of upper and lower surfaces of the beam are as follows u1(x,t)=−h[1− sin(πη0)]ψ(x,t) u2(x,t)= h[1+sin(πη0)]ψ(x,t) (2.2) Then, the longitudinal strain εx(x,y,t)= ∂u ∂x =−h { (η+η0) ∂2v ∂x2 − [sin(πη)+sin(πη0)] ∂ψ ∂x } (2.3) and the shear strain γxy(x,y,t)= ∂u ∂y + ∂v ∂x = πcos(πη)ψ(x,t) (2.4) The stresses – Hooke’s law σx(x,y,t) = E(y)εx(x,y,t) τxy(x,y,t)= G(y)γxy(x,y,t) (2.5) The simply supported beam with unsymmetrically varying mechanical properties of length L, depth h and width b is subjected to a uniformly distributed transverse load of intensity q or to axial compression force F0 (Fig. 3). 1166 K.Magnucki et al. Fig. 3. Scheme of the beam and loads The Hamilton principle t2∫ t1 [T − (Uε −W)] dt =0 (2.6) where: T is the kinetic energy, Uε – elastic strain energy, W – work of the load T = 1 2 bhρb L∫ 0 (∂v ∂t )2 dx Uε = 1 2 b L∫ 0 h/2∫ −h/2 [E(y)ε2x +G(y)γ 2 xy] dxdy W = L∫ 0 [ qv(x)+ 1 2 F0 (∂v ∂x )2] dx (2.7) and the equivalent – meanmass density of the beam ρb = 1 h h/2∫ −h/2 ρ(y) dy = 1 2 (ρ1+ρ2)= 1 2 ρ1(1+ √ e2) (2.8) Substitution of expressions (1.1)1 and (1.1)2 for the elasticity moduli and expressions (2.3) and (2.4) for strains into expression (2.7)2, after integration along depth of the beam, gives elastic strain energy as a functional of the two unknown functions Uε = 1 4 E1bh 3 L∫ 0 [ Cvv (∂2v ∂x2 )2 −2Cvψ ∂2v ∂x2 ∂ψ ∂x +Cψψ (∂ψ ∂x )2 +Cψ0 ψ2(x,t) h2 ] dx (2.9) where Cvv = 1 12 [ 1− 48 π2 η0+12η 2 0 + ( 1+ 48 π2 η0+12η 2 0 ) e2 ] Cψ0 = π2 4(1+ν1) (1+g2) Cvψ = 1 2π2 {(4−π2η0)[1− sin(πη0)]+(4+π2η0)[1+sin(πη0)]e2} Cψψ = 1 2 − sin(πη0)+sin2(πη0)+ [1 2 +sin(πη0)+sin 2(πη0) ] e2 Based on Hamilton’s principle (2.6) with consideration of expressions (2.7)1, (2.7)3 and (2.9), two differential equations of motion are obtained in the following form bhρb ∂2v ∂t2 + 1 2 E1bh 3 ( Cvv ∂4v ∂x4 −Cvψ ∂3ψ ∂x3 ) +F0 ∂2v ∂x2 = q Cvψ ∂3v ∂x3 −Cψψ ∂2ψ ∂x2 +Cψ0 ψ(x,t) h2 =0 (2.10) Bending, buckling and free vibration of a beam... 1167 The bendingmoment Mb(x)= b h/2∫ −h/2 yσx(x,y) dy (2.11) Substituting expression (2.5) for the normal stress, after integration along depth of the beam, one obtains the following equation Cvv d2v dx2 −Cvψ dψ dx =−2Mb(x) E1bh3 (2.12) It may be noticed that for static problems this equation is equivalent to equation (2.10)1. The position of the neutral axis is determined on the basis of the following condition – total axial force at the cross section h/2∫ −h/2 σx(x,y) dy =0 (2.13) Substituting expression (2.5) for the normal stress, after integration along depth of the beam, one obtains the following equation CNv d2v dx2 −CNψ dψ dx =0 (2.14) where CNv =(1+e2)η0− 2 π2 (1−e2) CNψ = 1 2 (1−e2)− (1+e2)sin(πη0) Based on this condition, the position of the neutral axis η0 = y0/h is obtained (Fig. 2). 3. Analytical solution of two differential equations of motion of the beam The system of two differential equations (2.10) for the beam is approximately solved with the use of two assumed functions v(x,t) = va(t)sin ( π x L ) ψ(x,t) = ψa(t)cos ( π x L ) (3.1) where: va(t), ψa(t) are functions of time t, which in the case of static problems become parame- ters. These functions satisfy the conditions of a simply supported beam. Substitution of functions (3.1) into equations (2.10) gives the following equations { bhρb d2va dt2 + 1 2 (π L )4 E1bh 3 [ Cvvva(t)− L π Cvψψa(t) ] − (π L )2 F0va(t) } sin ( π x L ) = q (π L )3 Cvψva(t)− (π L )2[ Cψψ + (λ π )2 Cψ0 ] ψa(t)= 0 (3.2) where λ = L/h is relative length of the beam. From equation (3.2)2, the function of time related to the shear effect is ψa(t)= π L kseva(t) (3.3) 1168 K.Magnucki et al. where the dimensionless coefficient of the shear effect is kse = Cvψ Cψψ + ( λ π )2 Cψ0 (3.4) It may be noticed that the value of this coefficient decreases with increasing relative length of the beam. Equation (3.2)1 with consideration of expression (3.3) is in the following form [ bhρb d2va dt2 + 1 2 (π L )4 E1bh 3(Cvv −kseCvψ)va(t)− (π L )2 F0va(t) ] sin ( π x L ) = q (3.5) and after application of Galerkin’s method is as follows bhρb d2va dt2 + 1 2 (π L )4 E1bh 3(Cvv −kseCvψ)va(t)− (π L )2 F0va(t)= 4 π q (3.6) This equation is the base for detailed studies of the bending, buckling and free vibration of the simply supported beamwith unsymmetrically varying mechanical properties. Condition (2.14) for calculation of the position of the neutral axis of the beam (Fig. 2) with consideration of functions (3.1) and (3.2)1 and expression (3.3) takes form of a transcendental equation η0−kse sin(πη0)− ( 2 π2 − 1 2 kse )1−e2 1+e2 =0 (3.7) It may be noticed that for large relative length of the beam (λ →∞, kse =0), the position of the neutral axis is determined as η (lim) 0 = η kse=0 0 = 2 π2 1−e2 1+e2 (3.8) Example. The following data of the beam are assumed: Poisson’s ratios: ν1 = ν2 = 0.33, e2 = 0.010, 0.025, 0.050, and relative length λ = 4,6, . . . ,14,∞. The dimensionless valu- es η0 (3.7) and η (lim) 0 (3.8) of the position of the neutral axis of the beam are specified in Table 1. Table 1.The dimensionless values η0 of the position of the neutral axis e2 λ 4 6 8 10 12 14 ∞ 0.010 0.2019 0.2001 0.1995 0.1992 0.1990 0.1989 0.1986 0.025 0.1962 0.1943 0.1936 0.1933 0.1931 0.1930 0.1928 0.050 0.1870 0.1850 0.1843 0.1839 0.1838 0.1836 0.1833 The graph of the dimensionless values η0 and η (lim) 0 of the position of the neutral axis of the beam is shown in Fig. 4. For the homogeneous beam (e2 =1), the neutral axis is located in the middle depth of the beam (η0 =0). Bending, buckling and free vibration of a beam... 1169 Fig. 4. The graph of dimensionless values η0 for the position of the neutral axis of the beam 4. Bending of the beam, static problem – analytical solution The simply supported beam with unsymmetrically varying mechanical properties is subjected to a uniformly distributed transverse load of intensity q (Fig. 3). On the basis of equation (3.6) for the static problem (d2va/dt 2 =0) and F0 =0, the relative maximum deflection is obtained ṽmax = vmax L = kv max qλ3 E1b (4.1) where the dimensionless coefficient of the maximal deflection is kv max = 8 π5(Cvv +kseCvψ) (4.2) In the case of large relative length of the beam (λ →∞, kse =0), this coefficient of themaximun deflection is k(lim)v max = 8 π5Cvv (4.3) Example. The following data of the beam are assumed: Poisson’s ratios: ν1 = ν2 = 0.33, e2 = 0.010,0.050, . . . ,0.50,1.0, and relative length λ = 5,10,15,20,∞. The values of the dimensionless coefficient of the maximum deflection kvmax and k (lim) v max are specified in Table 2. The graph of the values of the dimensionless coefficient of the maximum deflection kv max and k (lim) v max is shown in Fig. 5. For the homogeneous beam (e2 =1) of large relative length, the dimensionless coefficient of the maximum deflection k (lim) v max =48/π 5. 1170 K.Magnucki et al. Table 2.The dimensionless coefficient kvmax of the maximum deflection e2 λ 5 10 15 20 ∞ 0.01 0.6215 0.5978 0.5934 0.5919 0.5899 0.05 0.5314 0.5084 0.5042 0.5027 0.5008 0.10 0.4550 0.4329 0.4288 0.4274 0.4256 0.25 0.3312 0.3116 0.3080 0.3067 0.3051 0.50 0.2431 0.2267 0.2237 0.2226 0.2213 1.0 0.1733 0.1610 0.1587 0.1579 0.1569 Fig. 5. The graph of the values of the dimensionless coefficient of the maximum deflection 5. Buckling of the beam, static problem – analytical solution The simply supported beam with unsymmetrically varying mechanical properties is subjected to axial compressionwith the force F0 (Fig. 3). On the basis of equation (3.6) for static problem (d2va/dt 2 =0) and q =0, the critical force is obtained F0,CR = (π λ )2 kFCRE1bh (5.1) where the dimensionless coefficient of the critical force is kFCR = 1 2 (Cvv −kseCvψ) (5.2) In the case of large relative length of the beam (λ →∞, kse =0), this coefficient of the critical force is k (lim) FCR = 1 2 Cvv (5.3) Bending, buckling and free vibration of a beam... 1171 Example. The following data of the beam are assumed: Poisson’s ratios: ν1 = ν2 = 0.33, e2 =0.010,0.050, . . . ,0.50,1.0, and relative length λ =25,30,35,40,∞. The values of the dimensionless coefficient of the critical force kFCR and k (lim) FCR are specified in Table 3. The λ values in the buckling problem are larger than those for bending, since the critical loads for short beams would be very high and, therefore, elastic-plastic buckling would arise. Table 3.The values of the dimensionless coefficient kFCR and k (lim) FCR of the critical force e2 λ 25 30 35 40 ∞ 0.01 0.022112 0.022126 0.022135 0.022141 0.022159 0.05 0.026038 0.026058 0.026070 0.026077 0.026102 0.10 0.030629 0.030655 0.030671 0.030681 0.030714 0.25 0.042698 0.042742 0.042769 0.042787 0.042844 0.50 0.058845 0.058916 0.058959 0.058987 0.059078 1.0 0.082985 0.083091 0.083155 0.083197 0.083333 The graph of the values of the dimensionless coefficient of the critical force kFCR and k (lim) FCR is shown in Fig. 6. Fig. 6. The graph of the values of the dimensionless coefficient of the critical force For the homogeneous beam (e2 =1) of large relative length, the dimensionless coefficient of the critical force k (lim) FCR =1/12. 6. Free vibration of the beam, dynamic problem – analytical solution The simply supported beamwith unsymmetrically varying mechanical properties is not loaded (q =0, F0 =0) (Fig. 3). Equation (3.6) for the dynamic problem is as follows ρb d2va dt2 + 1 2 (π L )4 E1h 2(Cvv −kseCvψ)va(t)= 0 (6.1) 1172 K.Magnucki et al. The equation is solved with the use of the assumed function va(t)= va sin(ωt) (6.2) where: va is the amplitude of flexural vibration, ω – fundamental natural frequency. Substituting this function into equation (6.1), after simple transformation, one obtains the fundamental natural frequency ω = (π λ )2 kω √ E1 ρbh 2 (6.3) where the dimensionless coefficient of the fundamental natural frequency is kω = √ 1 2 (Cvv −kseCvψ)= √ kFCR (6.4) Taking into account expression (5.3), one formulates k (lim) ω = √ k (lim) FCR. Example. The following data of the beam are assumed: Poisson’s ratios: ν1 = ν2 = 0.33, e2 = 0.010,0.050, . . . ,0.80,1.0, and relative length λ = 5,10,15,25,∞. The values of the dimensionless coefficient of the fundamental natural frequency kω and k (lim) ω are specified in Table 4. Table 4.The values of the dimensionless coefficient kω and k (lim) ω of the natural frequency e2 λ 5 10 15 25 ∞ 0.01 0.14502 0.14787 0.14842 0.14870 0.14886 0.05 0.15683 0.16034 0.16101 0.16136 0.16156 0.10 0.16949 0.17376 0.17458 0.17501 0.17526 0.25 0.19866 0.20481 0.20601 0.20663 0.20699 0.50 0.23187 0.24011 0.24173 0.24258 0.24306 0.80 0.25984 0.26954 0.27146 0.27246 0.27303 1.0 0.27465 0.28497 0.28701 0.28807 0.28868 The graph of the values of the dimensionless coefficient of the fundamental natural frequency kω and k (lim) ω is shown in Fig. 7. In the case of the homogeneous beam (e2 = 1) of large relative length, the dimensionless coefficient of the critical force k (lim) FCR = √ 3/6. 7. Numerical calculations – FEM study 7.1. Numerical FEM model The numerical analysis of the beam with unsymmetrically varying mechanical properties is carried outwith thehelp of theSolidWorks software. The simulation assumed the samegeometry parameters andmechanical properties as those used in the analytical calculations. The beam is modelled using 3D finite elements in 20 layers, each with different mechanical properties satisfying expressions (1.1). Taking into account symmetry of the structure, a half of the beam is considered (Fig. 8). Bending, buckling and free vibration of a beam... 1173 Fig. 7. The graph of the values of the dimensionless coefficient of the fundamental natural frequency Fig. 8. Boundary conditions and load for the bending problem in the FEM study Therefore, the following boundary conditions are adopted: • for x =0 – the simple support – v(0) displacements in the y direction are zero; • for x = L/2 – themiddle of the beam – u(L/2) displacements in the x direction are zero. The numerical study of bending, buckling and free vibration is restrained to the xy-plane, similarly as in the case of the analytical approach. SolidWorks calculations have been carried out for beams with a rectangular cross-section of depth h =80mm, width b =20mm, and length values L = λh (400mm¬ L ¬ 3200 mm). 7.2. Bending of the beam, static problem – numerical FEM solution The beam is subjected to a uniformly distributed load of intensity q. A view to the bent half-beam is shown in Fig. 9. 1174 K.Magnucki et al. Fig. 9. Deflection of the beam (SolidWorks simulation) Results of the study are maximum deflections vmax [mm]. Based on expressions (4.1) and (4.2), values of the dimensionless coefficient kvmax are calculated. These values are specified in Table 5. Table 5. The values of the dimensionless coefficient kv max of the maximum deflection (FEM study) e2 λ 5 10 15 20 0.01 0.6136 0.5950 0.5886 0.5875 0.05 0.5248 0.5050 0.5007 0.5000 0.10 0.4496 0.4295 0.4257 0.4250 0.25 0.3272 0.3090 0.3062 0.3050 0.50 0.2400 0.2250 0.2222 0.2216 1.0 0.1696 0.1595 0.1580 0.1569 Values of the relative difference between analytical and FEM solutions are below 2.2%. The highest difference occurs for small relative length values λ. For greater λ values, the difference decreases. 7.3. Buckling of the beam, static problem – numerical FEM solution The half-beam under compression is shown in Fig. 10. Cross-sections and variation of me- chanical properties of the beam are the same as in the case of bending. Therefore, the buckling is analysed only in the xy-plane, similarly as for bending and free vibration. In order to avoid lateral buckling the z-displacements in the whole xy-plane are zeroed. Fig. 10. Boundary conditions and load for the buckling problem in the FEM study Buckling shape of the beam resulting from the numerical study is shown in Fig. 11. Results of the study are critical force values F0,CR [N]. Based on expressions (5.1) and (5.2), values of the dimensionless coefficient kFCR are calculated. These values are specified inTable 6. Values of the relative difference between analytical and FEM solutions are below 4.6%. Bending, buckling and free vibration of a beam... 1175 Fig. 11. Buckling shape of the beam (SolidWorks simulation) Table 6.The values of the dimensionless coefficient kFCR of the critical force (FEM study) e2 λ 25 30 35 40 0.01 0.0211 0.0213 0.0215 0.0216 0.05 0.0252 0.0255 0.0256 0.0258 0.10 0.0300 0.0302 0.0304 0.0306 0.25 0.0424 0.0429 0.0431 0.0435 0.50 0.0591 0.0598 0.0602 0.0607 1.0 0.0806 0.0811 0.0815 0.0817 7.4. Free vibration of the beam, dynamic problem – numerical FEM solution Free vibrations are computedwith the SolidWorks software for the FEMmodel composed of 10 layers of varyingmechanical properties.Ahalf-beam is adoptedwith the boundary conditions shown in Fig. 12. Fig. 12. Half-beammodel used in the free-vibration calculation Themiddle of the beam is placed in the left-hand side of the illustration. Hence, the x and z displacements are zeroed there. The right-hand part of the beam is simply supported and, therefore, the displacement y is blocked. The SolidWorks simulation tool used to compute the free-vibration frequencies provides angular frequencies ω of particular vibration modes. Nevertheless, in order to compare the analytical and numerical results, the dimensionless coefficient of the natural frequencies should be calculated for each case specified in Table 4. Taking into account expressions (2.8) and (6.3), one obtains 1176 K.Magnucki et al. kω = ωh (λ π )2 √ ρ1(1+ √ e2) 2E1 (7.1) In the case of the example presented in Fig. 13, the following data are assumed: e2 = 0.01, b = 20mm, h = 80mm, L = 400mm. Such a data set corresponds to the upper row and left- -hand column of Table 4. The angular frequency in this case is equal to ω =5093rad/s, which gives kω =0.14240. Fig. 13. An example result for e2 =0.01 and λ =5 The results kω of all the considered cases are presented in Table 7. Table 7.Values of the dimensionless coefficient kω of the natural frequency computed numeri- cally e2 λ 5 10 15 25 0.01 0.14240 0.14609 0.14680 0.14714 0.05 0.15530 0.15987 0.16075 0.16119 0.10 0.16775 0.17326 0.17434 0.17489 0.25 0.19658 0.20426 0.20578 0.20655 0.50 0.22962 0.23955 0.24153 0.24254 0.80 0.25743 0.26893 0.27125 0.27243 1.0 0.27214 0.28435 0.28681 0.28805 They perfectly comply with the results of Table 4 obtained analytically. The relative dif- ference values between analytical and FEM solutions are below 2%. The highest difference, equal to 1.8%, occurs for the example case mentioned above, whereas for the others they are significantly smaller. 8. Conclusions Theneutral axis of the studied beamdeviates from the geometric centre of the rectangular cross section as a result of unsymmetrical properties of the material. The location of the neutral axis is affected by the shear effect (Fig. 4). This effect is meaningful in the case of short beams and disappears for longer ones. The values of deflection, critical load and free-vibration frequencies depend on the position of the neutral axis. Therefore, in the case of the analytical approach, the Bending, buckling and free vibration of a beam... 1177 position of the axis should be determined first of all. The values of deflection, critical load and free-vibration frequencies obtained analytically have been compared to those computedwith the SolidWorks software. Itmay be noticed that the difference between both sets of the results does not exceed 5%. References 1. BirmanV., Byrd L.W., 2007,Modeling and analysis of functionally gradedmaterials and struc- tures,Applied Mechanics Reviews, 60, 195-216 2. Chen D., Kitipornchai S., Yang J., 2016a, Nonlinear free vibration of shear deformable san- dwich beamwith a functionally graded porous core,Thin-Walled Structures, 107, 39-48 3. Chen D., Yang J., Kitipornchai S., 2015, Elastic buckling, and static bending of shear defor- mable functionally graded porous beam,Composite Structures, 133, 54-61 4. Chen D., Yang J., Kitipornchai S., 2016b, Free and force vibrations of shear deformable functionally graded porous beams, International Journal of Mechanical Sciences, 108-109, 14-22 5. Debowski D., Magnucki K., Malinowski M., 2010, Dynamic stability of a metal foam rec- tangular plate, Steel and Composite Structures, 10, 2, 151-168 6. Feyzi M.R., Khorshidvand A.D., 2017, Axisymmetric post-buckling behaviour of saturated porous circular plates,Thin-Walled Structures, 112, 149-158 7. Jun L., Xiang H., Xiaobin L., 2016, Free vibration analyses of axially loaded laminated com- posite beams using a unified higher-order shear deformation theory and dynamic stiffness method, Composite Structures, 158, 308-322 8. Kapuria S., Bhattacharyya M., Kumar A.N., 2008, Bending and free vibration response of layered functionally gradedbeams:A theoreticalmodel and its experimental validation,Composite Structures, 82, 390-402 9. Kolakowski Z., Mania R.J., 2015, Dynamic response of thin FG plates with a static unsym- metrical stable postbuckling path,Thin-Walled Structures, 86, 10-17 10. Kubiak T., 2005, Dynamic buckling of thin-walled composite plates with varying width-wise material properties, International Journal of Solids and Structures, 45, 5555-5567 11. KubiakT., 2011,Estimation of dynamic buckling for composite columnswith open cross-sections, Computers and Structures, 89, 21-22, 2001-2009 12. Li L., Hu Y., Nonlinear bending and free vibration analyses of nonlocal strain gradient beams made of functionally gradedmaterial, International Journal of Mechanical Sciences, 107, 77-97 13. Magnucka-Blandzi E., 2011,Dynamic stability and static stress state of a sandwich beamwith ametal foam core using three modified Timoshenko hypothesis,Mechanics of Advanced Materials and Structures, 18, 2, 147-158 14. Mahi A., Bedia E.A.A., Tounsi A., 2015, A new hyperbolic deformation theory for bending and free vibration analysis of isotropic, functionally graded, sandwich and laminated composite plates,Applied Mathematical Modelling, 39, 2489-2508 15. Mojahedin A., Jabbari M., Khorshidvand A., Eslami M., 2016, Buckling analysis of func- tionally graded circular plates made of saturated porous material based on higher order shear deformation theory,Thin-Walled Structures, 99, 83-90 16. Sayyad A.S., Ghugal Y.M., 2017, Bending, buckling and free vibration of laminated composite and sandwich beams: A critical review of literature,Composite Structures, 171, 486-504 17. Song M., Kitipornchai S., Yang J., 2017, Free and force vibrations of functionally graded polymer composite plates reinforced with graphene nanoplatelets, Composite Structures, 159, 579-588 1178 K.Magnucki et al. 18. SmyczynskiM.J.,Magnucka-Blandzi E., 2018,Three-point bending of a sandwichbeamwith two binding layers – Comparison of two nonlinear hypotheses,Composite Structures, 183, 96-102 19. ThaiH.-T.,VoT.P., 2012,Bendingand free vibrationof functionally gradedbeamsusingvarious higher-order shear deformation beam theories, International Journal of Mechanical Sciences, 62, 57-66 20. Thai H.-T., Vo T.P., 2013, A new sinusoidal shear deformation theory for bending, buckling, and vibration of functionally graded plates,Applied Mathematical Modelling, 37, 3269-3281 21. Zenkour A.M., 2006, Generalized shear deformation theory for bending analysis of functionally graded plates,Applied Mathematical Modelling, 30, 67-84 22. ZhangY.,Wang S., Loughlan J., 2006, Free vibration analysis of rectangular composite lami- nates using a layerwise cubic B-spline finite strip method,Thin-Walled Structures, 44, 601-622 Manuscript received September 21, 2017; accepted for print June 7, 2018