Quasi-3D Hyperbolic Shear Deformation Theory for the Free Vibration Study of Honeycomb Microplates with Graphene Nanoplatelets-Reinforced Epoxy Skins

A novel quasi-3D hyperbolic shear deformation theory (QHSDT) with five unknowns is here employed, together with the Hamilton’s principle and the modified couple stress theory (MCST) to analyze the vibrational behavior of rectangular micro-scale sandwich plates resting on a visco-Pasternak foundation. The sandwich structure features a Nomex or Glass phenolic honeycomb core, and two composite face sheets reinforced with graphene nanoplatelets (GPLs). The effective properties of both face sheets are evaluated by means of the Halpin-Tsai and extended rule of mixture (ERM) micromechanical schemes. The governing equations of the problem are derived by applying the Hamilton’s principle, whose solutions are determined theoretically according to a classical Navier-type procedure. A parametric study checks for the effect of different material properties, length-scale parameters, foundation parameters and geometrical properties of the honeycomb cells, and the reinforcing GPLs, on the vibration response of the layered structure, which can be of great interest for many modern engineering applications and their optimization design.


Introduction
In the last decades, lightweight mechanical components and layered structures have increased the attention of many researchers and scientists, due to the increased demand in modern engineering, together with a possible reduction in their production cost. Among them, sandwich structures can be regarded as subset of multilayered composite structures consisting of outer facings and a soft core in-between, including foam, honeycomb, corrugated core, various bio-inspired cores, etc. The choice of sandwich materials depends on the structural functionality as well as on the lifetime loading, availability and cost. For example, Graphite-epoxy and carbon-epoxy multilayered facings are typically used in aerospace applications, whereas glass-epoxy or glass-vinyl ester are adopted in civil and marine layered structures. At the same time, the core of sandwich aerospace structures is often made of aluminum or many works from the literature have applied CPTs, FSDTs, or higher-order-theories (HSDTs) for the study of plate and shell structures even with complicated materials and geometries. For example, Khoa et al. [17] applied a HSDT to examine the vibration response of FG carbon nanotubes (CNTs) reinforced composites cylindrical shells in thermal environment. The same problem was also studied by Ibrahim et al. [18], according to FSDT, and coupled with thermal conditions. Li et al. [19] used CPT to model clamped honeycomb sandwich panels to study the nonlinear forced vibrational response. Many further applications of the HSDT to coupled problems of sandwich panels and shell structures can be found in [20][21][22][23][24][25][26][27][28]. A valid theoretical alternative to handle the plate structures is represented by the quasi-3D hyperbolic shear deformation theory (QHSDT) which accounts for both transverse shear and normal deformations and satisfies the zero traction boundary conditions on the plate surfaces without using any shear correction factor. In QHSDT the number of unknown functions involved in displacement field is only equal to five, instead of six or more unknowns required by the other shear and normal deformation theories. The computational efficiency of this method was recently verified in Refs. [29][30][31]. Inspired by these few pioneering works from the literature, in the present paper we propose a QHSDT to study the free vibration response of sandwich structures with a honeycomb core resting on a visco-Pasternak foundation. The governing equations of the problem are derived from the Hamilton's principle and solved in closed form via the Navier's method. The analytical solutions from our formulation are verified with those reported in literature, where a parametric investigations aims at determining the effect of the material variation, GPLs gradient index and dispersion patterns, geometry, internal cells angle, or thickness of layers on the natural frequencies for the selected sandwich structure.

Theoretical Formulation
Consider a rectangular sandwich plate with thickness h, length a, width b, as illustrated in Figure 1, together with the reference coordinate system (x, y, z). The sandwich structure is immersed within a visco-Pasternak elastic foundation, and it is made of a honeycomb core with thickness h c and two composite face sheets with thickness h t and h b at the top and bottom side, respectively. This means that the total thickness of the structure is h = h c + h t + h b . In the last relation, ϕ is an additional displacement variable that accounts for the effect of normal stress; g(z) and f(z) are expressed by the following functions [29] ( ) ( ) where f ′ (z) denotes the first derivative of function f with respect to z. The strain-displacement relations follow the von-Karman's assumptions [31]  In the current study, a QHSDT is adopted to define the position of an arbitrary point in the micro-model. The major advantage of using such a displacement field is that the problem is not limited to plane-strain conditions (i.e., ε zz 0), as typically occurs in the other 2-D theories such as FSDT, that could cause possible discrepancies between the theoretical and experimental results. Based on a QHSDT, the displacement field is defined as [32] U(x, y, where u and v stand for the displacement components along the x and y directions, respectively; w s , w b and w st are the transverse displacement components due to bending, shear and stretching effects, respectively, with w st (x, y, z, t) = g(z)ϕ(x, y, t) In the last relation, ϕ is an additional displacement variable that accounts for the effect of normal stress; g(z) and f (z) are expressed by the following functions [29] where f '(z) denotes the first derivative of function f with respect to z. The strain-displacement relations follow the von-Karman's assumptions [31] whereas the constitutive equations for the honeycomb core and FG-GPLs face sheets, read as follows [33] where C ij are the elastic constants for each part of the sandwich structure. More specifically, for the honeycomb core, the elastic constants read as follows [29] C 11c = where δ = 1 − 2ν 12 ν 13 ν 32 − ν 12 ν 21 − ν 13 ν 31 − ν 23 ν 32 (8) and Molecules 2020, 25, 5085 5 of 21 In the relations above, the Young's modulus, shear modulus, density, and Poisson's ratio are defined in a homogenized form by means of the mechanical properties E h , G h , ρ h and v h of the honeycomb material [34]. Besides, φ 0 is the internal cells angle of the honeycomb structure; ϕ 0 = h 0 /l 0 and γ 0 = t 0 /l 0 stand for the internal aspect ratio and dimensionless cells thickness, respectively, in which h 0 , l 0 and t 0 are the geometrical parameters defining the hexagonal cells as represented in Figure 1. For FG-GPLs reinforced face sheets, the elastic constants C ijf are given by The mechanical properties for both face sheets vary throughout the thickness, and they are clearly function of the effective material properties, defined, in turn, by means of the Halpin-Tsai micromechanical model, as follows [35] E(z) = 3 8 In the last relation, E M denotes the Young's modulus of the matrix; V GPL refers to the volume fraction of GPLs; ζ L , ζ W , η L and η W are the geometrical properties of GPLs, i.e., L GPL , h GPL , w GPL and E GPL being the length, thickness, width and Young's modulus of GPLs, respectively. It is noteworthy that the summation of GPLs and matrix volume fractions equals one, where the GPLs volume fraction is determined as where ρ GPL and ρ M refer to the density of the reinforcement phase and matrix, respectively. Moreover, g GPL is the weight fraction of GPLs that obey the following relations for three different dispersion patterns though the face sheets thicknesses [36] For a parabolic pattern For a linear pattern (24) in which the positive and negative signs are related to the top and bottom face sheets, respectively. For a uniform pattern In Equations (23)-(25), λ P , λ L and λ U are the gradient index of GPLs for their parabolic, linear, and uniform dispersion patterns, referred to the total GPLs content, as reported in Table 1. Table 1. Graphene nanoplatelets (GPLs) gradient index for different values of their total content [36].

Total GPLs Content (Percentage
The further properties for the face sheets are the Poisson's ratio and density, which are determined via the ERM as [37] ρ

Governing Equations of the Problem
The Hamilton's principle is here applied to gain the governing equations of the problem [38] where Π, Λ and K denote the applied external work, the strain energy, and the kinetic energy for the sandwich structure, respectively. The strain energy of the system consists of two parts: the classical strain energy and the energy component from the MCST. The following relation is used to define the total strain energy for the selected sandwich structure [39] where m ij and χ ij stand for the higher-order stresses and symmetric rotation gradient tensor, respectively, defined in the following where l m is the MCST material length scale parameter, and µ is the Lame's parameter. Moreover, the components of the symmetric rotation gradient tensor can be determined using the following compact relation This means that in which, the infinitesimal rotation vector Θ is defined as which means In addition, the kinetic energy for the whole microstructure can be defined as [40].
where U, V and W refer to the displacement components introduced in Equation (1). For a structure resting on a visco-Pasternak elastic foundation, the external work due to the substrate can be defined as follows [41] where K W is the Winkler parameter, K G is the shear layer parameter, and C d denotes the damping parameter, respectively. By substitution of Equations (29), (35), (36) into the Hamilton's principle (28), after a mathematical manipulation we get the following governing equations of motion in terms of displacement field δu : Molecules 2020, 25, 5085 8 of 21 δv : −C 120 Molecules 2020, 25, 5085 9 of 21 δw s : δϕ : More details about the coefficients in Equations (37)- (41), are reported in the Appendix A.

Analytical Solution Procedure
The differential equations of the Equations (37)-(41) are solved analytically according to the Navier's procedure in this section. Therefore, for a simply supported structure, we consider the following theoretical expressions for the displacement components [42] u(x, y, t) = U cos(αx) sin(βy)e iωt , v(x, y, t) = V sin(αx) cos(βy)e iωt , w b (x, y, t) = W b sin(αx) sin(βy)e iωt , w s (x, y, t) = W s sin(αx) sin(βy)e iωt , ϕ(x, y, t) = Φ sin(αx) sin(βy)e iωt (42) in which U, V, W s , W b and Φ are the unknown coefficients. In addition, α and β are defined as mπ/a and nπ/b, respectively, where m and n are the mode numbers along the length and width direction, respectively. After substituting Equation (42) into Equations (37)-(41), the equations of motion gain the following compact form

Numerical Results
In this section we illustrate the numerical results, in terms of vibration response, for a microsandwich plate with a honeycomb core made of Nomex or Glass phenolic, and Epoxy-reinforced GPLs as face sheets. The Nomex has the following properties: E s = 3.2 GPa, ρ = 48 kg/m 3 , and ν = 0.4. For the Glass phenolic, the same properties of Ref. [43] are assumed herein. The Epoxy matrix and GPLs reinforcement phase for the face sheets have the following properties [44] E GPL = 1.01 TPa, ρ GPL = 1062.5 kg/m 3 , ν GPL = 0.186, L GPL = 2.5 µm, The microplate has a total height equal to 150 µm, 80% of whose total height corresponds to the core, and the rest is equally divided between the two face sheets. The length of the square plate is ten-fold of its thickness. The internal cell angle, aspect ratio, and dimensionless cells thickness are assumed to be π/6, 1, and 0.1, respectively. Moreover, the material length-scale parameter is kept as 15 µm according to Ref. [36].
To check for the reliability of our formulation, we compare the results for a single-layer FG-GPL-reinforced square microplate with predictions by Thai et al. [45]. The comparative results are summarized in Table 2, in terms of dimensionless natural frequencies defined as Ω = ωa 2 /h ρ M /E M , for various mode numbers, while considering the effect of the aspect ratio (a/h) and length-scale parameter-to-total thickness (l m /h). Based on Table 2, a very good agreement is observable between the results from our formulation and those ones from Ref. [45], where some negligible differences are related to the different kinematic assumptions, and/or different solution techniques.
In Table 3, we also summarize the natural frequencies for different mode numbers, as computed according to a MCST or a classical elasticity theory (CET), for a varying internal cell angle from 30 • up to 60 • . Based on results in Table 3, note that the mode number and internal cell angle yield a reverse effect on the natural frequency of the sandwich microplate, whereby an increasing mode number and a decreasing internal cell angle get higher values of the natural frequencies. It seems also that CET-based predictions are always more conservative than those once based on a MCST, in agreement with findings in Refs. [46][47][48][49][50] from the literature. Another key aspect of the problem can be the sensitivity of the response to various GPLs dispersions in the Epoxy matrix over a wide range of mode numbers, as listed in Table 4. It is worth noticing that the sandwich structure becomes stiffer for an increased quantity of GPLs as reinforcing phase, and the natural frequency enhances dramatically in each mode number. It seems also that a linear dispersion of GPLs in the Epoxy matrix with λ L = 2 is more effective than other types of distribution with λ P = λ U = 1 for an overall increase in the structural stiffness. This shows that the GPLs dispersion coefficient plays a crucial role, more than the type of GPLs dispersion, for an increase in the natural frequency. Figure 2 shows the variation of the natural frequency for the sandwich microplate against the l m /h ratio, for different dispersions of GPLs. By increasing l m /h rational value, and keeping constant the total thickness of the sandwich model, the natural frequency increases monotonically, for each fixed value of λ L , λ P , λ U . This behavior is due to a reduced flexibility of the sandwich microplate which corresponds to a stiffness and stability enhancement. For each type of GPLs dispersion, a higher distribution coefficient obtains higher natural frequencies.  Figure 3 plots the effect of the aspect ratio, a/b, on the natural frequency of the microstructure for different GPLs dispersion coefficients. For each fixed GPLs dispersion coefficient and type, an increased aspect ratio up to one clearly reduces the natural frequency reaching the minimum value for a cubic sandwich structure. Once this minimum value is passed, the aspect ratio rolling up causes a monotonic increase in the natural frequency for each selected GPLs dispersion coefficient and type. A further systematic analysis is also performed to check for the sensitivity of the natural frequency alternation with the lm/h ratio, under different assumptions for the honeycomb core material in Figure  4. Based the plots in this figure, it is worth observing that the most rigid sandwich microstructure is obtained for a uniform Nomex honeycomb core material, where the most flexible one is reached for an Epoxy/Glass Phenolic core material. All the other results based on an Epoxy/Nomex honeycomb or Uniform/Glass Phenolic core material assumption are very close to each other, and fall always within the previous two cases. As also plotted in Figure 5, the natural frequency decreases monotonically for an increasing geometrical ratio hGPL/LGPL of the reinforcing phase, as predicted by a CET or a MCST, respectively, while assuming three different rational values for LGPL/WGPL, namely, LGPL/WGPL = 1; 5/3; 2. This means that the GPLs length variations (reduction or enhancement) have a direct relationship with the natural frequency, stiffness and rigidity. For each selected theory, an increased value of LGPL/WGPL reduces gradually the natural frequency for each fixed value of hGPL/LGPL. Based on a comparative evaluation of the curves in Figure 5, it can be noted that MCST provides always a higher natural frequency compared to the CET.  Figure 3 plots the effect of the aspect ratio, a/b, on the natural frequency of the microstructure for different GPLs dispersion coefficients. For each fixed GPLs dispersion coefficient and type, an increased aspect ratio up to one clearly reduces the natural frequency reaching the minimum value for a cubic sandwich structure. Once this minimum value is passed, the aspect ratio rolling up causes a monotonic increase in the natural frequency for each selected GPLs dispersion coefficient and type. A further systematic analysis is also performed to check for the sensitivity of the natural frequency alternation with the l m /h ratio, under different assumptions for the honeycomb core material in Figure 4. Based the plots in this figure, it is worth observing that the most rigid sandwich microstructure is obtained for a uniform Nomex honeycomb core material, where the most flexible one is reached for an Epoxy/Glass Phenolic core material. All the other results based on an Epoxy/Nomex honeycomb or Uniform/Glass Phenolic core material assumption are very close to each other, and fall always within the previous two cases. As also plotted in Figure 5, the natural frequency decreases monotonically for an increasing geometrical ratio h GPL /L GPL of the reinforcing phase, as predicted by a CET or a MCST, respectively, while assuming three different rational values for L GPL /W GPL , namely, L GPL /W GPL = 1; 5/3; 2. This means that the GPLs length variations (reduction or enhancement) have a direct relationship with the natural frequency, stiffness and rigidity. For each selected theory, an increased value of L GPL /W GPL reduces gradually the natural frequency for each fixed value of h GPL /L GPL . Based on a comparative evaluation of the curves in Figure 5, it can be noted that MCST provides always a higher natural frequency compared to the CET. Molecules 2020, 25, x FOR PEER REVIEW 5 of 25     In Figure 6 we analyze the effect of the viscoelastic foundation on the vibration response, while providing the 3D plot of the natural frequency for different combinations of KW, KG under three different assumptions for the damping parameter Cd = 500; 1000; 1500 (N.s/m) is provided. By increasing the Winkler and Pasternak parameters (KW, KG) the structural stiffness increases together with the natural frequency for each fixed value of Cd. Based on the three plots, it is worth mentioning the great damping effect on the frequency response, where a decreased value of Cd obtains higher frequencies for each fixed combination of KW, KG. In Figure 6 we analyze the effect of the viscoelastic foundation on the vibration response, while providing the 3D plot of the natural frequency for different combinations of K W , K G under three different assumptions for the damping parameter C d = 500; 1000; 1500 (N·s/m) is provided. By increasing the Winkler and Pasternak parameters (K W , K G ) the structural stiffness increases together with the natural frequency for each fixed value of C d . Based on the three plots, it is worth mentioning the great damping effect on the frequency response, where a decreased value of C d obtains higher frequencies for each fixed combination of K W , K G .  The effect of the internal aspect ratio φ 0 and dimensionless cell thickness γ 0 on the first natural frequency of the sandwich microplate is plotted in Figure 7. Based on the results in this figure, larger magnitudes of γ 0 lead to an increased system stability. On the other hand, a clear reduction in the structural stiffness and frequency is gained by internal aspect ratio enhancement and honeycomb core thickness reduction in the case of fixed internal cells angle equal to 30 • . This means that, for a constant value of total thickness, a lower face sheet thickness to core thickness ratio results in a higher stiffness and weaker flexibility. structural stiffness and frequency is gained by internal aspect ratio enhancement and honeycomb core thickness reduction in the case of fixed internal cells angle equal to 30°. This means that, for a constant value of total thickness, a lower face sheet thickness to core thickness ratio results in a higher stiffness and weaker flexibility.
Moreover, based on the curvatures plotted in Figures 8 and 9 which represent the first natural frequency versus the honeycomb core internal cell angle θ0 for different internal aspect ratios φ0, it seems that an enhancement of both parameters gets a natural frequency reduction. In addition, Figure  9 illustrates that the thicker honeycomb core provides higher structural stiffness and natural frequency. As a final parametric investigation, we check for the variation of the first natural frequency with the lm/h , based on the MCST or CET, under the assumption of three different core thicknesses. Based on the plots in Figure 10, it should be noted that the natural frequency increases significantly for higher values of lm/h ratio, when the problem is tackled by a MCST, whereas it remains almost unaffected by lm/h according to a CET. This confirms, once again, the great importance of adopting a size-dependent approach instead of classical formulations. Moreover, based on the curvatures plotted in Figures 8 and 9 which represent the first natural frequency versus the honeycomb core internal cell angle φ 0 for different internal aspect ratios ϕ 0 , it seems that an enhancement of both parameters gets a natural frequency reduction. In addition, Figure 9 illustrates that the thicker honeycomb core provides higher structural stiffness and natural frequency. As a final parametric investigation, we check for the variation of the first natural frequency with the l m /h, based on the MCST or CET, under the assumption of three different core thicknesses. Molecules 2020, 25, x FOR PEER REVIEW 8 of 25

Figure 9.
Effect of the thickness ratio on the first natural frequency of the structure. Figure 9. Effect of the thickness ratio on the first natural frequency of the structure.
Based on the plots in Figure 10, it should be noted that the natural frequency increases significantly for higher values of l m /h ratio, when the problem is tackled by a MCST, whereas it remains almost unaffected by l m /h according to a CET. This confirms, once again, the great importance of adopting a size-dependent approach instead of classical formulations.

Conclusions
In this work, a QHSDT is employed to investigate the vibrational behavior of sandwich honeycomb microplates with two GPLs' composite face sheets, resting on elastic foundations. The equations of motion are obtained by applying the Hamilton's principle, where the Navier-type solutions are determined in analytical form. Based on a large systematic investigation, it is noted that an Epoxy/Nomex honeycomb core makes the sandwich structure less flexible than Epoxy/Glass phenolic and uniform glass phenolic core materials, whereby a uniform Nomex honeycomb core provides the highest structural stiffness. Moreover, a larger dimensionless cell thickness (γ0) yields an increased stability in the system, whereas internal aspect ratio elevation provides structural stability reduction along with the system's stiffness and natural frequency. The results based on a MCST are compared to predictions from CET to provide a clear understanding about vibrational responses' sensitivity to size-dependent parameters. In agreement with findings from the literature, a CET always produces more conservative results compared to an MCST, which justifies the necessity of adopting non classical approaches instead of the classical ones. The proposed model together with our numerical results could be very useful for the design and manufacturing of many aerospace, automotive or shipbuilding engineering applications, where honeycomb structures are recommended for their great capability to tolerate high pressures and stresses despite their light structure.

Conclusions
In this work, a QHSDT is employed to investigate the vibrational behavior of sandwich honeycomb microplates with two GPLs' composite face sheets, resting on elastic foundations. The equations of motion are obtained by applying the Hamilton's principle, where the Navier-type solutions are determined in analytical form. Based on a large systematic investigation, it is noted that an Epoxy/Nomex honeycomb core makes the sandwich structure less flexible than Epoxy/Glass phenolic and uniform glass phenolic core materials, whereby a uniform Nomex honeycomb core provides the highest structural stiffness. Moreover, a larger dimensionless cell thickness (γ 0 ) yields an increased stability in the system, whereas internal aspect ratio elevation provides structural stability reduction along with the system's stiffness and natural frequency. The results based on a MCST are compared to predictions from CET to provide a clear understanding about vibrational responses' sensitivity to size-dependent parameters. In agreement with findings from the literature, a CET always produces more conservative results compared to an MCST, which justifies the necessity of adopting non classical approaches instead of the classical ones. The proposed model together with our numerical results could be very useful for the design and manufacturing of many aerospace, automotive or shipbuilding engineering applications, where honeycomb structures are recommended for their great capability to tolerate high pressures and stresses despite their light structure.

Acknowledgments:
The authors would like to thank Ehsan Arshid for his tireless efforts and guidance during this research.

Conflicts of Interest:
The authors declare no conflict of interest.