An Alternative Electro-Mechanical Finite Formulation for Functionally Graded Graphene-Reinforced Composite Beams with Macro-Fiber Composite Actuator

With its extraordinary physical properties, graphene is regarded as one of the most attractive reinforcements to enhance the mechanical characteristics of composite materials. However, the existing models in the literature might meet severe challenges in the interlaminar-stress prediction of thick, functionally graded, graphene-reinforced-composite (FG-GRC)-laminated beams that have been integrated with piezoelectric macro-fiber-composite (MFC) actuators under electro-mechanical loadings. If the transverse shear deformations cannot be accurately described, then the mechanical performance of the FG-GRC-laminated beams with MFC actuators will be significantly impacted by the electro-mechanical coupling effect and the sudden change of the material characteristics at the interfaces. Therefore, a new electro-mechanical coupled-beam model with only four independent displacement variables is proposed in this paper. Employing the Hu–Washizu (HW) variational principle, the precision of the transverse shear stresses in regard to the electro-mechanical coupling effect can be improved. Moreover, the second-order derivatives of the in-plane displacement parameters have been removed from the transverse-shear-stress components, which can greatly simplify the finite-element implementation. Thus, based on the proposed electro-mechanical coupled model, a simple C0-type finite-element formulation is developed for the interlaminar shear-stress analysis of thick FG-GRC-laminated beams with MFC actuators. The 3D elasticity solutions and the results obtained from other models are used to assess the performance of the proposed finite-element formulation. Additionally, comprehensive parametric studies are performed on the influences of the graphene volume fraction, distribution pattern, electro-mechanical loading, boundary conditions, lamination scheme and geometrical parameters of the beams on the deformations and stresses of the FG-GRC-laminated beams with MFC actuators.


Introduction
Composite structures with sensors and actuators, which are called smart structures, are used in a wide range of engineering applications [1][2][3]. Among the various available smart materials, the piezoelectric material is one of the most commonly used materials that can be used as both actuators and sensors. Macro-fiber composites (MFCs) are a new type of piezoelectric intelligent materials that were developed by NASA. The rectangular cross-section of the piezoelectric fibers in the MFC increases the contact area between the interdigitated electrodes and piezoceramic fibers in order to improve the transfer of electricity [4][5][6]. With a higher electro-mechanical conversion efficiency and better flexibility, the MFC has broad application prospects and research value.
In order to achieve the vibration control and shape simulation of smart structures with MFC patches, Zhang et al. [7] studied the active shape and vibration control of Static, buckling and dynamic analyses of laminated composite plates with MFC actuators have been extensively investigated [8][9][10][11][12][13][50][51][52][53][54]. In addition, piezoelectric materials were used in the smart structures for shape, vibration and buckling control [7,12,[55][56][57]. However, the research on the structural behaviors of FG-GRC-laminated plates with piezoelectric actuators is quite limited. Mao and Zhang [58,59] analyzed the linear and nonlinear vibration and buckling behaviors of piezoelectric FG-GRC plates. Using the transformed-differential-quadrature method, Malekzadeha et al. [60] studied the free vibration of eccentric annular FG-GRC plates that were embedded in piezoelectric layers. Nguyena et al. [61] proposed an isogeometric finite-element formulation for the bending and transient analysis of porous FG-GRC plates that were embedded in piezoelectric layers. The free vibration of FG-GRC sandwich plates that were enclosed by piezoelectric layers was investigated by Majidi-Mozafari et al. [62]. The wave-characteristic analysis of cylindrical FG-GRC sandwich structures with piezoelectric surface layers was reported by Li and Han [63,64]. Khayat at al. [65] investigated the uncertainty propagation in the nonlinear dynamic responses of smart, porous cylindrical FG-GRC sandwich shells. Guo et al. [66] investigated the nonlinear dynamic behavior of three-phase composite plates that were made with MFCs in the polymer and with graphene skins.
Owing to the electro-mechanical coupling effect and the discontinuity of the material properties throughout the thickness of the FG-GRC-laminated structures with piezoelectric layers, transverse shear deformations play an important role. In order to precisely predict the transverse shear stresses of FG-GRC-laminated structures with MFC actuators, the compatible conditions of transverse shear stresses at the interfaces must be fulfilled. As can be observed in the above literature, FSDT and HSDT combined with different solution methods are used for the mechanical analysis of FG-GRC structures with piezoelectric layers. However, FSDT and HSDT are unable to satisfy the compatible conditions of transverse shear stresses at the interfaces. As a result, FSDT and HSDT may fail to model the transverse shear stresses of the FG-GRC-laminated plates with MFC actuators. Therefore, there is a genuine requirement to propose an efficient and accurate model with which to study this issue.
From the above literature survey, it was found that the comprehensive interlaminar shear-stress analysis of the thick FG-GRC-laminated beams with MFC actuators are scarce in the open literature. This paper aims to fill this research gap by developing an efficient finite-element formulation based on an attractive electro-mechanically coupled model. The proposed model fulfills the compatible conditions of transverse shear stresses at the interfaces. By using the 3D elasticity equations and the HW variational principle [67,68], the accuracy of the transverse shear stresses in terms of the electro-mechanical coupling effect is improved. Due to the second-order derivatives of the in-plane displacement, variables are removed from the interlaminar shear-stress components, such that the finite-element formulation can be easily developed. Thus, based on the proposed electro-mechanical model, a simple C 0 -type two-node beam element is developed for the interlaminar-stress analysis of thick FG-GRC-laminated beams with MFC actuators. The performance of this finite-element formulation is appraised by comparing it with the 3D elasticity solutions and the results are calculated from the chosen models. Comprehensive parametric studies are conducted in order to explore the influences of distribution pattern, volume fraction, lamination sequence, as well as geometric parameters of the beams on the deformations and stresses of FG-GRC-laminated beams with MFC actuators. The research findings will provide essential information for the reliable design of smart FG-GRC-laminated structures with great potential for engineering applications.

Effective Material Properties of GRC
An FG-GRC-laminated beam integrated with an MFC layer is shown in Figure 1. The length of the beam is a. The thicknesses of the FG-GRC-laminated beam and the MFCactuator layer are h and h p , respectively. Each GRC layer may have a different graphene volume fraction. The graphene reinforcement distribution is functionally graded and of the piece-wise type in terms of thickness. The four distribution patterns of the FG-GRC beam are shown in Figure 2, denoted as UD, FG-V, FG-O and FG-X. Different colors are used to represent the distributions of the volume fraction in terms of thickness. A dark color indicates a large volume fraction, whereas a light color indicates a small volume fraction. First of all, the graphene volume fractions of the FG-V GRC beams are set to [(0.11) 2 /(0.09) 2 /(0.07) 2 /(0.05) 2 /(0.03) 2 ]. In this type of FG-GRC beam, the volume fraction of the top layer is the largest and the volume fraction of the bottom layer is the smallest. Secondly, the graphene volume fractions are arranged as [0.11/0.09/0.07/0.05/0.03] s for a ten-layer beam and expressed as FG-X. Different from the FG-V type, the top and bottom GRC layers have the largest graphene volume fraction and the layers near the mid-surface have the smallest graphene volume fraction. Then, the arrangement of graphene volume fractions [0.03/0.05/0.07/0.09/0.11] s is referred to FG-O. In this type, the volume fraction is maximum at the mid-surface and minimum at the top and bottom layers. Finally, for comparative purposes, the type of UD GRC beams with the same thickness and graphene volume fraction V G = 0.07 are considered.  Based on the extended Halpin-Tsai model [38], the effective elastic moduli of the GRC layer are given by in which the superscripts G and m denote graphene and matrix, a G , b G and h G are the effective length, width and thickness of the graphene sheet, respectively. The volume fractions V G and V m satisfy the condition V G + V m = 1. The graphene efficiency parameters η i (i = 1, 2, 3) in Equation (1) can be obtained by the MD simulation. The density and Poisson ratio of the composite media can be expressed as where v G 12 and v m 12 denote the Poisson ratio of graphene sheet and matrix, respectively.

Displacement Field
The displacement field of the present beam theory is a combination of the first-order shear-deformation theory and the efficient zigzag functions, and is given by where superscript k represents the kth ply of the laminated beam, u k and w k denote the assumptions of in-plane displacement and transverse displacement for the kth ply of the beam, u 0 is the in-plane axial displacement, u 1 is the bending rotation, u z is the amplitude function and w 0 is the transverse displacement. The zigzag function Φ k (z) in Equation (4) can be explicitly written as where u k are the interfacial axial displacements.
, and z k is defined in Figure 1.
In terms of the displacement field given by Equation (4), the transverse shear strain is expressed as where The transverse shear stress for the kth ply is expressed as (8) in which Q k 44 represent the shear modulus of the kth ply In Equation (8), it is found that the transverse shear stress contains two parts: the continuity conditions are applied only on the zigzag-dependent part, τ k xz , that is, on the transverse shear stresses obtained by removing the shear measures η x . The continuity conditions are expressed as Equation (9) provides N-1 conditions for the interfacial axial displacements. The other two conditions can be obtained from the zero-value conditions on the outer plate surfaces, which are given by where t and b are used to denote the bottom and top layers of the beam. According to the Equations (9) and (10), the u k given by Equation (5) are obtained. Based on the assumed displacement field, the strains are shown as where The constitutive relations for the kth ply of a cross-ply FG-GRC-laminated beam are given by where Q k ij are the transformed elastic stiffness coefficients of the kth ply. For the piezoelectric layer, the stress components are given by [69] where e k ij are the transformed piezoelectric modulus of the kth ply. The electric field E can be obtained by the following equation where ξ denotes the electrostatic potential, which is assumed as in which φ(x) represents the axial distribution of the applied electrostatic potential, V t and h p denote the applied electric voltage and the thickness of the MFC layer, respectively.

Improved Transverse Shear Stresses Considering Electro-Mechanical Coupling Effect
Owning to the electro-mechanical coupling effect and the sudden change of material characteristics at the interfaces of the FG-GRC-laminated beams with piezoelectric layers, the transverse shear stresses that are calculated directly from the constitutive equations have low accuracy. To achieve an improved through-the-thickness variation of transverse shear stresses in terms of the electro-mechanical coupling effect, the 3D equilibrium equations and HW variational principle are used. The key feature of this formulation is that the assumed transverse shear stresses have a superior accuracy.
For the laminated beam, the 3D equilibrium equation can be reduced to the 2D equilibrium equation. Disregarding the influence of the body force, the 2D equilibrium equation is expressed as Integrating Equation (18) with respect to the z-coordinate and enforcing the condition τ k xz (−h/2) = 0 yields the following form for the transverse shear stress Substituting Equations (11)- (17) into Equation (19), yields where By applying the condition τ k xz (h/2 + h p ) = 0, the second-order derivative function ∂ 2 u 0 /∂x 2 is removed from the Equation (20), yielding where in which The assumed transverse shear stress τ k xz in Equation (23) can a priori fulfill the interlaminar continuity conditions at the interfaces and the zero-traction conditions on the outer beam surfaces. The transverse shear stresses contain the second-order derivatives of the in-plane displacement variables, which will lead to difficulties in finite-element modeling. Thus, the HW variational principle [67,68] is employed to eliminate the secondorder derivatives of the in-plane displacement parameters. In terms of the HW variational principle, a function for the laminated beams can be stated as where u, ε and σ are the independent variables that are subjected to variation, and Ω is the region occupied by a whole beam.
Π ext (u) = S qudS (29) in which q is a traction vector, and the transverse shear strain γ k xz corresponding to the assumed transverse shear stress τ k xz is given by Evaluation of the first variation of the Equation (27) and equating to zero results in Applying integration by parts, the following equation is given by where D and D S are the differential operators. Introducing Equation (33) into Equation (32) yields According to the second integral in Equation (34), the following equation can be given by Introducing Equations (11), (23) and (30) into Equation (35), results in Using the above equation, the second-order derivatives of the displacement parameters U x can be rewritten as where Introducing Equation (37) into Equation (23) yields the final expression of τ k

Finite-Element Formulation
In this section, a two-node beam element is used in order to study the FG-GRClaminated beams with MFC actuators. Accordingly, the displacement variables for each element can be expressed as where N 1 = (1 − ξ)/2, N 2 = (1 + ξ)/2; ξ = (x − x c )/l e , x c = x i + l e /2, l e = x i+1 − x i ; x i and x i+1 are defined in Figure 3. By means of Equations (11)-(13), the strain vector is written as T and d e i are the degrees of freedom at the ith node. The degrees of freedom at each node are given by According to the Equation (33), the element-stiffness matrix of the proposed beam element can be obtained. Equation (33) within the ith element can be rewritten as where W i represents the work of external forces. Introducing Equations (11)-(17), (39) into Equation (45), the element-stiffness matrix is expressed as (46) in which Q is the material-constant matrix, the strain matrix B is given in Equation (44), and the strain matrixB is expressed aŝ where R k xz = 1/Q k 44 , M k j are given in Equation (39), and j = 1~3. By means of the following equation, the vector of nodal displacement d e can be expressed as where F e and F e p are the loading vector and the applied-electric-potential vector for an element, respectively.
in which Assembling all the element-stiffness matrices, the global equation is obtained in the following form: where K is the stiffness matrix; F and F p are the global loading vector and the appliedelectric-potential vector, respectively.

Numerical Results and Discussion
In this section, various numerical examples are presented in order to assess the proposed finite-element formulation. Additionally, the influences of the electro-mechanical loading, graphene distribution pattern, volume fraction, lamination scheme, thickness of the MFC-actuator layer, and length-to-thickness ratio on the static behaviors of FG-GRClaminated beams with MFC actuators are carried out.

Verification Study
As a typical example, a simply supported laminated beam (0 • /90 • /0 • ) subjected to a sinusoidal mechanical load q = q 0 sin(πx/a) is firstly analyzed. The material characteristics of the laminated beam are [70]: 25. Dimensionless results can be obtained by the relations: Convergence studies are first conducted, and the results with different mesh densities are presented in Table 1. It is found that the converged results agree well with the 3D elasticity solutions (Exact) [70]. In addition, the displacements and stresses converge at the mesh density of 48, which may be considered to be adequate for the thick laminated beam with a piezoelectric layer. Table 2 presents the displacements and stresses obtained from different models. The present results are compared with the exact solutions [70] and the results obtained from Reddy's higher-order shear-deformation theory (HSDT-R), trigonometric higher-order shear-deformation theory (HSDT-S&G), first-order sheardeformation theory (FSDT) and classical lamination theory (CLT). The results of HSDT-R, HSDT-S&G, FSDT and CLT are provided by Sayyad and Ghugal [71]. In Table 2, it can be observed that the present model can produce more accurate results when compared with other models.  In order to further assess the present finite-element formulation, the simply supported FG-GRC-laminated beam with an MFC layer (0 • /90 that was subjected to a sinusoidal load q = q 0 sin(πx/a) is investigated. The material properties for the matrix and graphene are [72]: 177. The graphene efficiency parameters [33] for different volume fractions are shown in Table 3. The assumptions G 23 = G 13 = 0.5G 12 are made for the numerical analysis [33]. The material constants for the MFC piezoelectric layer [50] Figures 4-7 were computed by the authors in terms of Pagano's approach [70]. It can be seen from Figures 4-7 that the present plate theory can accurately yield the stresses without using any post-processing procedure. However, Reddy's higher-order theory (HSDT-R) [73], four-variable refined higher-order theories (HSDT-T [42], HSDT-A [40]) and first-order shear-deformation theory (FSDT) are apparently less accurate in predicting the transverse shear stresses. It is also observed that the HSDT-R, HSDT-T, HSDT-A and FSDT models are unable to satisfy the compatible conditions of transverse shear stresses. Therefore, for the interlaminar-stress analysis of FG-GRC-laminated beams with MFC layers, the compatible conditions of transverse shear stresses at the interfaces should be satisfied.

Parametric Study
In this section, parametric studies on the static analysis of FG-GRC-laminated beams with MFC actuators under electro-mechanical loads are conducted. The influences of the graphene volume fraction, distribution pattern, electro-mechanical loading, thickness of the MFC-actuator layer, lamination scheme, and the length-to-thickness ratio on the deformations and stresses of an FG-GRC-laminated beam with an MFC actuator under combined sinusoidal, mechanical and electrical loads are presented. Unless otherwise specified, the beams are simply supported, and the thicknesses of the FG-GRC-laminated beam and the MFC layer are h = 10 mm and hp = 0.5 mm, respectively. The lamination scheme is taken to be (0 • /90 • /0 • /90 • /0 • /0 • /90 • /0 • /90 • /0 • /M). The material characteristics of the matrix, graphene and MFC piezoelectric layer are the same as those used in the previous section. Dimensionless results can be obtained by the following relations: Firstly, the FG-GRC-laminated beam with an MFC actuator under electro-mechanical loading (q 0 = 100, V = 100) is investigated. Figure 8 depicts the normalized transverse deflections of the UD GRC plates with MFC-actuator layers. It is found that the transverse deflections are greatly influenced by the graphene volume fraction V G . The transverse deflections can be significantly reduced by adding a very small amount of graphene nanofillers. Therefore, graphene can help to improve the stiffness of the beam and suppress the deflections. Table 4 presents the transverse displacement and stresses of thick, laminated FG-GRC beams with MFC layers (a/h = 8) for the different graphene distributions and boundary conditions under a lateral load (q 0 = 100). The following three boundary conditions (BCs) are considered: two ends of the beam are simply supported (SS), one end (x = 0) of the beam is clamped and the other end is simply supported (CS), and two ends of the beam are clamped (CC). In Table 4, it is found that the transverse displacement and stresses are strongly influenced by the boundary conditions. The more constraint that is applied to the ends of the beam, the more rigidity of the boundaries that is obtained. The results in Table 4 indicate that when other parameters are the same, the SS beam has the highest deflection and the CC beam has the lowest owing to the clamped end having more rigidity than the simply supported end.  The through-the-thickness variations of the in-plane displacements and stresses of the thick FG-GRC-laminated beam (a/h = 4) with an MFC actuator under electro-mechanical loading (q 0 = 100, V = 100) for different graphene distribution patterns are shown in Figure 9. It can be seen from Figure 9 that the graphene distribution pattern has significant effects on the in-plane displacements and stresses of smart FG-GRC-laminated beams. Different graphene distribution patterns lead to different stiffnesses of the beam, which then lead to different displacements and stresses. Consequently, the deformations and stresses of smart FG-GRC-laminated beams can be regulated by changing the internal distribution pattern in engineering applications. In addition, it is also observed from Figure 9 that when a voltage is applied, the transverse shear stresses at the interface between the MFC piezoelectric layer and the FG-GRC beams will suddenly increase.
The normalized displacements and stresses of simply supported FG-GRC-laminated beams with MFC actuators under a sinusoidal mechanical load (q 0 = 100) and electric loads (V = 0, V = 100 and V = −100) are shown in Table 5. The span-to-thickness ratio a/h is taken to be 10. In Table 5, it is found that the reversal of deflections and stresses can be observed due to the change in voltage polarity. In addition, the results in Table 5 reveal that the actuation responses of the displacements and stresses are very large. The effects of the thickness of the MFC-actuator layer on the actuation responses of a smart, thick FG-GRC-laminated beam (a/h = 5) subjected to an electro-mechanical loading (q 0 = 100, V = 100) are shown in Table 6. The results in Table 6 indicate that the actuation effects on the deformations and stresses are decreased with the increase in the thickness of the MFC-actuator layer h p . Based on the above analysis, it can be concluded that the electro-mechanical loading and the activated MFC piezoelectric layer have remarkable impacts on the static analysis of FG-GRC-laminated beams with MFC actuators.   Table 7 presents the effect of the length-to-thickness ratio (a/h) on the deflection and transverse shear stress of FG-X and FG-V GRC beams with an MFC actuator subjected to an electro-mechanical loading (q 0 = 100, V = 100). It is found that the actuating response of the thick FG-GRC-laminated beams with MFC actuators is more significant than that of the thin beams. Additionally, it should be noted that the τ xz (h/2) is the transverse shear stress between the MFC piezoelectric layer and the FG-GRC beam. It can be seen from Table 7 that the actuation responses of the transverse shear stresses are decreased with the increase in a/h.

Conclusions
Based on an attractive electro-mechanical coupled-beam model, a C 0 -type finiteelement formulation was developed for the interlaminar shear stress analysis of a thick FG-GRC-laminated beam with an MFC actuator. Using the HW variational principle, the improved transverse shear stresses in terms of the electro-mechanical coupling effect were derived. In order to assess the accuracy and efficiency of the present formulation, the results acquired from other models and the 3D elasticity solutions were also utilized for comparison. Furthermore, comprehensive parametric studies were carried out in order to analyze the effects of various parameters on the deformations and stresses. By analyzing typical static problems, some conclusions can be drawn: (1) The present model can satisfy the compatible conditions of transverse shear stresses, and the results acquired from the present model agreed well with the 3D elasticity solutions.
(2) The stiffness of the FG-GRC-laminated beam with an MFC actuator can be significantly enhanced by adding a small amount of graphene. The displacements and stresses of the FG-GRC-laminated beam with an MFC layer are sensitive to the graphene distribution pattern. Therefore, the transverse shear stresses can be optimized by adjusting the graphene distribution pattern.
(3) The actuating responses of the thick FG-GRC-laminated beams with MFC actuators are more significant than that of the thin beams. When the electric voltages were applied, the maximum transverse shear stress was obtained at the interface between the MFC piezoelectric layer and the GRC laminates.
In addition, the effects of interlaminar stresses on the stability and dynamic behaviors of smart FG-GRC-laminated structures will be conducted in future work.  Institutional Review Board Statement: Not applicable.

Conflicts of Interest:
The authors declared no potential conflict of interest with respect to the research, authorship, and/or publication of this article.