**1. Introduction**

Composites are made up of two or more constituent materials that, when combined, generate a material with physical and chemical characteristics different from the individual components. Composite structures have presented outstanding advances in the last decades. The recent applications of fiber-reinforced composite structures to marine [1], civil [2], aerospace [3] and other industries [4–6] have already been reviewed in the open literature. Fast progress in manufacturing has resulted to the necessity for the improvement of composites in terms of strength, stiffness, density, and lower cost.

Graphene is a carbon allotrope, the thinnest known material, i.e., one carbon atom thick, which additionally is exceptionally strong—almost 200 times stronger than steel [7–9]. Furthermore, graphene presents outstanding electric [10,11], heat [12] and optical [13] properties. The excellent properties of graphene make it an excellent candidate as reinforcing material in composites. Consequently, many studies can be found in the literature investigating experimentally [14–17], computationally [18–20], or analytically [21,22] the effect of graphene in nanocomposites. The superb behavior of graphene nanocomposites has driven the research and technology for finding advanced applications. Chang and

Wu [23] have reported recent energy-related development of graphene nanocomposites in solar energy conversion, such as photovoltaic and photoelectrochemical devices and artificial photosynthesis, as well as electrochemical energy devices and improvements in environmental applications of functionalized graphene nanocomposites for the detection and removal of heavy metal ions, organic pollutants, gas and bacteria. Kumar et al. [24] have systematically examined the recent progress in illustrating the complexities of mechanical and thermal behavior of graphene- and modified graphene-based polymer nanocomposites. They also reported the potential applications, existing challenges, and potential perspectives regarding the multi-scale capabilities and promising advancements of the graphene family-based nanocomposites materials.

The fiber and matrix-dominant behavior of fiber-based polymer composites are crucial in many mechanical, aerospace, or structural applications. Pre-mixing the polymer matrix with nanoparticles could enhance the through-thickness or matrix-dominant properties. Rafiee et al. [25] employed several graphene-based nanomaterials, including graphene oxide and its thermally reduced version (rGO), graphene nanoplatelets, and multi-walled carbon nanotubes to modify the epoxy matrix and the surface of glass fibers. Unmodified and modified epoxy and fibers were used for fabricating multiscale glass fiber-reinforced composites, following the vacuum-assisted resin transfer molding process. The composites obtained combined improvements in both the fiber and matrix-dominant properties, resulting in superior composites. Kostagiannakopoulou et al. [26] developed a new class of carbon fiber-reinforced polymers with a nano-modified matrix, based on graphene nano-species. Their results have shown that nano-doped composites present a significant improvement of the interlaminar strain energy release rate of the order of 50% for the case of graphene nano-platelets. Ta¸s and Soykok [27] have determined engineering constants of carbon nanotube based unidirectional carbon fiber-reinforced composite lamina theoretically with two different approaches. They have modeled using the finite element method and analyzed the behavior of a specific stacking sequence and configuration of a composite plate under static loadings. Their results have shown a considerable improvement of the plate stiffness due to the presence of nanotubes.

Although some works have already been done on vibrations of composite structures reinforced by nanoparticles, there are only a few works, which focused on vibration behavior of carbon fiber-based laminate composites with pure graphene inclusions. Hulun et al. [28] have investigated free vibration of graphene nanoplatelet (GPL) reinforced laminated composite quadrilateral plates using the element-free IMLS-Ritz method. They concluded that increasing GPLs weight fraction can increase the natural frequency of quadrilateral plate. Shen et al. have investigated nonlinear vibration of functionally graded graphene-reinforced composite laminated plates [29], functionally graded graphene-reinforced composite laminated cylindrical panels resting on elastic foundations [30], and functionally graded graphene-reinforced composite laminated cylindrical shells [31] in thermal environments. Wang et al. [32] have proposed a two-dimensional elasticity model of laminated graphene-reinforced composite (GRC) beams, where the graphene disperses uniformly in each layer, but the graphene volume fraction may vary from layer to layer. It is found that the laminated GRC beam with graphene distribution pattern X has the least deflection and highest fundamental frequency at extreme length-to-thickness ratios, but it has the highest deflection and least fundamental frequency at a very low length-to-thickness ratio due to its reduced transverse shear stiffness.

In the present paper, a computational approach is proposed for the estimation of vibration characteristics of carbon fiber laminate composite plates with graphene inclusions. The method firstly computes the size-dependent mechanical properties of graphene used in the composite. In the second phase, the orthotopic properties of the hybrid matrix are calculated adopting Halpin-Tsai models for specific volume fractions of graphene in the matrix. Knowing the elastic constants of hybrid matrix, the elastic mechanical properties of the composite lamina can be evaluated. In the next step, finite element models are developed for various configurations of laminated composite plates using conventional shell elements. Applying different boundary conditions and solving the free vibration problem, modes of vibration and corresponding natural frequencies are finally extracted. The effect of

graphene volume fraction on vibration characteristics of various composite plates is determined. To the authors' best knowledge, it is the first time where this holistic and relatively simple approach (four-stage technique), taking into account the size-dependent behavior of graphene inclusions, is proposed for the vibration analysis of graphene-carbon fiber-reinforced composite structures. *Materials* **2020**, *13*, x FOR PEER REVIEW 3 of 15 (four-stage technique), taking into account the size-dependent behavior of graphene inclusions, is proposed for the vibration analysis of graphene-carbon fiber-reinforced composite structures.

#### **2. Computational Approach 2. Computational Approach**

#### *2.1. Problem Definition 2.1. Problem Definition*

Figure 1 depicts the composite plate considered to be analyzed. In Figure 1a, the finite element model of the graphene in nanoscale is illustrated. The information about the size-dependent mechanical performance of graphene derived from nanoscale is used for the determination of the elastic constants of hybrid matrix (Figure 1b) using Halpin-Tsai homogenization model. The elastic constants of the hybrid matrix are utilized in order to compute the orthotropic mechanical properties of the unidirectional composite lamina (Figure 1c). Then, the finite element modeling of a composite plate considering a specific stacking sequence (Figure 1d) can be performed for its vibration analysis. Figure 1 depicts the composite plate considered to be analyzed. In Figure 1a, the finite element model of the graphene in nanoscale is illustrated. The information about the size-dependent mechanical performance of graphene derived from nanoscale is used for the determination of the elastic constants of hybrid matrix (Figure 1b) using Halpin-Tsai homogenization model. The elastic constants of the hybrid matrix are utilized in order to compute the orthotropic mechanical properties of the unidirectional composite lamina (Figure 1c). Then, the finite element modeling of a composite plate considering a specific stacking sequence (Figure 1d) can be performed for its vibration analysis.

#### *2.2. Graphene Mechanical Properties in Nanoscale 2.2. Graphene Mechanical Properties in Nanoscale*

function [33]

The elastic mechanical behaviour of different-sized graphene can be numerically examined and predicted using a spring-based finite element approach [7,33]. According to those approaches, threedimensional, two-nodded, spring-based, finite elements of three degrees of freedom per node are combined in order to model effectively the interatomic interactions presented inside the graphene. The calculated variations of mechanical elastic constants have been approached by appropriate sizedependent non-linear functions of two independent variables, i.e., length and width for graphene, in order to express the analytical rules governing the elastic behaviour. The numerical results have been already validated through comparisons with corresponding data given in the open literature. The elastic mechanical behaviour of different-sized graphene can be numerically examined and predicted using a spring-based finite element approach [7,33]. According to those approaches, three-dimensional, two-nodded, spring-based, finite elements of three degrees of freedom per node are combined in order to model effectively the interatomic interactions presented inside the graphene. The calculated variations of mechanical elastic constants have been approached by appropriate size-dependent non-linear functions of two independent variables, i.e., length and width for graphene, in order to express the analytical rules governing the elastic behaviour. The numerical results have been already validated through comparisons with corresponding data given in the open literature.

The Young's modulus in TPa of the graphene can be predicted using the following rational

(1) *nm* 10 ݈ ,ݓ൏,0

ଵାభ௪ାభାమ௪మାమమାమ௪

ܧ൫ݓ, ݈൯ = బାబభ௪ାబభାబమమାబమ௪

The Young's modulus in TPa of the graphene can be predicted using the following rational function [33]

$$E\_{\mathcal{S}^{\mathcal{T}}}(w\_{\mathcal{S}^{\mathcal{T}}}l\_{\mathcal{S}^{\mathcal{T}}}) = \frac{P\_0 + a\_{01}w\_{\mathcal{S}^{\mathcal{T}}} + b\_{01}l\_{\mathcal{S}^{\mathcal{T}}} + b\_{02}l\_{\mathcal{S}^{\mathcal{T}}} + c\_{02}w\_{\mathcal{S}^{\mathcal{T}}\mathcal{S}^{\mathcal{T}}}}{1 + a\_1w\_{\mathcal{S}^{\mathcal{T}}} + b\_1l\_{\mathcal{S}^{\mathcal{T}}} + a\_2w\_{\mathcal{S}^{\mathcal{T}}} + b\_2l\_{\mathcal{S}^{\mathcal{T}}} + c\_2w\_{\mathcal{S}^{\mathcal{T}}\mathcal{S}^{\mathcal{T}}}}, 0 < w\_{\mathcal{S}^{\mathcal{T}}}l\_{\mathcal{S}^{\mathcal{T}}} \le 10 \text{ nm} \tag{1}$$

where *P*0, *a*01, *b*01, *b*02, *c*02, *a*1, *b*1, *a*2, *b*2, and *c*<sup>2</sup> are constants, which are determined by non-linear fitting, while *wgr* and *lgr* are the graphene width and length, respectively. The values of the constants of Equation (1), concerning both the two chirality directions of graphene are presented in Table 1. The accuracy of this equation compared to the finite element predictions and expressed by the coefficient of determination is calculated to be over than 99.5%.

**Table 1.** Fitting parameters of Equation (1) for the prediction of Young's modulus variation in TPa.


#### *2.3. Graphene*/*Polymer Matrix Physical Properties*

The three-dimensional (3D) Halpin-Tsai model for randomly oriented discontinuous rectangular reinforcements is utilized to compute the Young's modulus of the hybrid graphene/polymer matrix. According to this model, the Young's modulus of the graphene/polymer matrix can be predicted by the following equation [34]

$$E\_{m-gr} = \frac{1}{5}E\_L + \frac{4}{5}E\_T \tag{2}$$

where

$$E\_L = \frac{1 + \xi\_L \eta\_L V\_{\mathcal{S}^r}}{1 - \eta\_L V\_{\mathcal{S}^r}} E\_m \tag{3}$$

$$E\_T = \frac{1 + \xi\_T \eta\_T V\_{\mathcal{S}^r}}{1 - \eta\_T V\_{\mathcal{S}^r}} E\_m \tag{4}$$

and

$$\eta\_L = \frac{\frac{E\_{\text{ST}}}{E\_m} - 1}{\frac{E\_{\text{ST}}}{E\_m} + \xi\_L}, \eta\_T = \frac{\frac{E\_{\text{ST}}}{E\_m} - 1}{\frac{E\_{\text{ST}}}{E\_m} + \xi\_T} \tag{5}$$

$$\xi\_L = 2\frac{l\_{\mathcal{gr}}}{t\_{\mathcal{gr}}} + 40V\_{\mathcal{gr}}{}^{10}, \quad \xi\_T = 2\frac{w\_{\mathcal{gr}}}{t\_{\mathcal{gr}}} + 40V\_{\mathcal{gr}}{}^{10} \tag{6}$$

where *Vgr* and *tgr* stand for volume fraction and thickness of graphene, respectively. Furthermore, *E<sup>m</sup>* is the Young's modulus of the polymer matrix, while ξ*T*, ξ*L*, η*L*, and η*<sup>T</sup>* are Halpin-Tsai parameters, depending on the geometry of the reinforcement.

Assuming that Poisson's ratio *v<sup>m</sup>* for both pure polymer matrix and hybrid matrix are approximately similar, it is calculated that Poisson's ratio of graphene/polymer matrix is

$$
\upsilon\_{m-gr} = \upsilon\_m \tag{7}
$$

Consequently, the shear modulus of the hybrid matrix exhibiting quasi-isotropic behaviour may be evaluated from the next equation

$$\mathcal{G}\_{m-\mathcal{g}^r} = \frac{E\_{m-\mathcal{g}^r}}{2(v\_{m-\mathcal{g}^r} + 1)}\tag{8}$$

Considering the inertia effect for the dynamic problem, the equivalent mass density of the hybrid matrix can be determined using rule of mixtures as

$$
\rho\_{m-gr} = V\_{gr}\rho\_{gr} + V\_m\rho\_m \tag{9}
$$

where ρ*gr*, ρ*m*, and *V<sup>m</sup>* are the mass density of graphene, mass density of matrix, and the volume fraction of the matrix, respectively. It is also known that

$$V\_{\mathcal{R}^r} + V\_m = 1\tag{10}$$

#### *2.4. Physical Properties of Composite Lamina*

Considering a unidirectional lamina, its Young's modulus in longitudinal direction, E1, and Poisson's ratio, υ<sup>12</sup> can be computed by employing the standard rule of mixture [35] as follows:

$$E\_1 = E\_{f11} \ V\_f + E\_{m-gr} \ V\_{m-gr} \tag{11}$$

$$
\omega\_{12} = \upsilon\_f V\_f + \upsilon\_{m-gr} V\_{m-gr} \tag{12}
$$

where *Ef*11, *V<sup>f</sup>* , υ*<sup>f</sup>* and *Vm-gr* are Young's modulus of fiber along in longitudinal direction, volume fraction of fiber, Poisson's ratio of fiber and volume fraction of the hybrid matrix, respectively. The remaining orthotropic material properties, i.e., Young's modulus in transverse direction, *E*2, Poisson's ratio, υ23, and shear moduli, *G*<sup>12</sup> and *G*23, can be computed by the next general Halpin-Tsai [36] expression

$$\frac{p}{P\_{m-gr}} = \frac{1 + \xi \eta V\_f}{1 - \eta V\_f} \tag{13}$$

where *Pm-gr* are the related properties of hybrid matrix. Moreover, η is an experimental factor, which is computed by using the next equation,

$$\eta = \frac{\frac{P\_f}{P\_{m-gr}} - 1}{\frac{P\_f}{P\_{m-gr}} + \xi} \tag{14}$$

where *P<sup>f</sup>* are the related properties of fiber. Notice that ξ is a measure of reinforcement geometry, which affected by loading conditions. Here, ξ is chosen as 2 for computation of *E*<sup>2</sup> and as 1 for computations of υ23, *G*<sup>12</sup> and *G*23.

Considering inertia effects, the equivalent mass density of the composite structure can be determined using rule of mixtures as

$$
\rho\_c = V\_f \rho\_f + V\_{m-gr} \rho\_{m-gr} \tag{15}
$$

where ρ*<sup>f</sup>* , ρ*m*−*gr* and *Vm*−*gr* are the mass density of the fiber, mass density of the hybrid matrix, and volume fraction of the hybrid matrix, respectively. It is also known that

$$V\_f + V\_{m-gr} = 1\tag{16}$$
