Nonlinear low-velocity impact response of GRC beam with geometric imperfection under thermo-electro-mechanical loads

This paper presents an investigation on thermo-electro-mechanical nonlinear low-velocity impact behaviors of the geometrically imperfect functionally graded (FG) graphene-reinforced composite (GRC) beam with surface-bonded piezoelectric layers. Both uniformly distributed and FG patterns of graphene nanoplatelets (GPLs) are considered along the thickness of the GRC host beam. The effective Young’s modulus is calculated by the Halpin–Tsai model. The Poisson ratio and mass density are calculated by the rule of mixture. The modified nonlinear Hertz contact law is employed to predict the impact force between the spherical impactor and the geometrically imperfect GRC piezoelectric beam during impacting. By considering the first-order shear deformation theory and von Kármán nonlinear displacement–strain relationship, the nonlinear governing equations are obtained by Hamilton principle and dispersed by differential quadrature (DQ) method. The Newmark-β method associated with Newton–Raphson iterative process is adopted to parametrically identify the impact force and the dynamic response of the system. The effects of geometric imperfection, weight fraction and distribution pattern of GPLs, temperature variation, thickness of piezoelectric layer and impactor’s initial velocity on nonlinear low-velocity impact behaviors of geometrically imperfect GRC beams are discussed in detail. Our results illustrate that the coupling effect of geometric imperfection and thermo-electro-mechanical load has a significant effect on the nonlinear low-velocity impact behavior of GRC beam, and GPLs distributing into the piezoelectric layers is better for reducing the impact response of geometrically imperfect GRC piezoelectric beam.

structures [1][2][3], especially for the vibration, buckling and nonlinear dynamic behaviors [4][5][6]. Sheinman et al. [6] studied the imperfection sensitivity of an isotropic beam under a nonlinear elastic foundation. Farokhi et al. [7] investigated the effect of geometric imperfection on the nonlinear resonance of the Timoshenko microbeam. Wadee [8] proposed the one-dimensional imperfection model to simulate various global and localized imperfections existing in engineering structure, and then conducted a detailed research on buckling and post-buckling behaviors of the geometrically imperfect sandwich panel. Their results shown that the mechanical behaviors of structures are very sensitive to the geometric imperfections. Obviously, the research on the dynamic behaviors of structures with geometric imperfections is of great significance in engineering design.
Except the isotropic homogeneous structures mentioned above, more and more lightweight and high strength materials and high-performance composite structures are emerging in the field of aerospace and engineering applications, such as the functionally graded materials (FGM), fiber-reinforced composite structures and carbon-based composite structures [9][10][11][12]. In the framework of isogeometric analysis [13], Nguyen and his co-workers calculated the static and dynamic behaviors of the high-performance composite structures, such as the vibration of FGM piezoelectric porous microplates [14] and postbuckling of functionally graded carbon nanotube-reinforced composite (CNTRC) shells [15]. When the high-performance composite structures operate under thermo-electro-mechanical loads, their mechanical behaviors can significantly be affected by their working conditions. Yang and his co-authors [16,17] analyzed the large amplitude vibration and thermal bifurcation buckling of carbon nanotubereinforced composite (CNTRC) beams bonded with piezoelectric layers. Liew's group [18,19] researched active vibration control of FG composite plates bonded with the piezoelectric actuator and sensor. Wu et al. [20,21] analyzed the post-buckling and free vibration of FG-CNTRC beam with geometric imperfections under thermal-mechanical-electrical loads.
Recently, graphene nanoplatelets (GPLs) have shown tremendous potentials in the development of high-performance composites due to its unique mechanical properties since its first report by Novoselov et al. [22]. Extensive theoretical and experimental investigations have shown that the dispersion of GPLs in the polymer matrix can significantly improve its mechanical properties [23][24][25]. Yang's group firstly proposed the functionally graded (FG) GPLs-reinforced composite (GRC) structures and carried out a detailed study on their bending, buckling and vibration behaviors [26][27][28][29]. Mao et al. [30][31][32][33] conducted a detailed analysis of the static and dynamic behaviors of FG-GRC structures. Rahimi et al. [34,35] employed a semi-analytical solution to obtain the three-dimensional bending and free vibration behaviors of FG-GRC cylindrical shells. Furthermore, many scholars have studied the vibrations and stability of FG-GRC structures under various loads, including transverse excitation [36,37], axial compression [38], thermal loads [39] and impact loads [40]. Results illustrated that the GPL reinforcements can significantly enhance the stiffness of the GRC structures, and the enhancement effect largely depends on the distribution pattern of GPLs in matrix.
According to the researches about nonlinear behaviors of FG-GRC structures with geometric imperfections, shapes and amplitudes of geometric imperfections have significant effects on the buckling [41,42], nonlinear vibration [4,43] and resonance [44] of FG-GRC structures. It is worth noting that impact is one of common loads for engineering structures, which may result in internal damage and even lead to severe failure of structures. Introducing the modified Hertz model, Fan et al. [45,46] discussed the lowvelocity impact response of FG-GRC structures resting on visco-elastic foundations. Based on the energybalance and spring-mass model, Dong et al. [47] predicted the low-velocity impact response of FG-GR cylindrical shells under axial load and thermal loads. Selim et al. [48] examined the impact behaviors of FG-GRC plates resting on Winkler-Pasternak elastic foundations by using a meshless approach.
The existing literature have researched the lowvelocity impact behaviors of intact GRC structures and the nonlinear dynamic behaviors of geometrically imperfect GRC structures in detail separately. To our best knowledge, there is no relevant reports on the impact behavior of GRC beam with geometric imperfections. But it is the utmost important for the safety and longevity of engineering structures, especially for the high-performance composite structures working under complex loading, such as the situation among thermal, electrical and mechanical loads.
In order to investigate the couple effect of the lowvelocity impact and geometric imperfections on the nonlinear dynamic behaviors of the GRC beam subjected to thermal-mechanical-electrical loads, the present work introduces a theoretical model and numerical analysis method. For the geometric imperfect GRC beam with surface-bonded piezoelectric layers subjected to low-velocity impact and thermoelectro-mechanical loads, the governing equations of motion are derived based on the first-order shear deformation theory, von Kármán nonlinear displacement-strain relationship and Hamilton's principle, and dispersed by differential quadrature (DQ) method. The Newmark-b method associated with Newton-Raphson iterative method is employed to study the effects of the weight fraction and distribution pattern of GPLs, the types and amplitude of geometric imperfections, temperature variation, impact velocity, as well as thickness of the piezoelectric layer on impact force and impact response of the geometrically imperfect GRC piezoelectric beam.
2 Theoretical models 2.1 Geometry model Figure 1a illustrates a GRC beam with geometric imperfection. The length, width and thickness of the GRC beam are L, b and h, respectively. Two piezoelectric layers with a thickness of h p are perfectly bonded on the upper and lower surfaces of the geometrically imperfect GRC beam and respectively applied with voltage U 0 . The hinged supported boundary conditions are considered. The GRC beam is subjected a uniform temperature variation DT ¼ T À T 0 and it is stress-free at initial temperature T 0 . It is assumed that the material properties of GRC beam are independent of the varying temperatures to focus on the coupled effect of the temperature variation and geometric imperfection on the nonlinear low-velocity impact response of the GRC piezoelectric beam. The geometrically imperfect GRC piezoelectric beam is impacted by a spherical impactor with initial velocity V 0 . m i and R i respectively denote the mass and radius of the impactor.
The geometric imperfection w Ã x ð Þ is simulated in the form of the product of the trigonometric and hyperbolic functions [8] where a is a constant, representing the localization degree of the imperfection, and A 0 is the dimensionless amplitude of the geometric imperfection. The spherical impactor impacts to the concave surface of the GRC beam for A 0 [ 0, as shown in Fig. 1, and the convex surface is impacted for A 0 \ 0. Note that w Ã ðxÞ is symmetric about x=L ¼ c, and b is the half-wave numbers of the imperfection in x-axis. Sine, global (G) and localized (L) imperfections are considered. As presented in Fig

Nonlinear Hertz model
Based on Newton's second law, the governing equation of motion for the impactor can be expressed as where y i t ð Þ is the displacement of the impactor, and F C t ð Þ is the impact force between the impactor and the impacted geometrically imperfect GRC piezoelectric beam.
When the contact area between the impactor and the impacted beam is much smaller than the geometries of the beam, Hertz contact theory can be utilized to calculate the impact force. The nonlinear lowvelocity impact process can be divided into two phases, namely, loading phase and unloading phase [46,49]. The loading phase starts from the moment that the impactor impacts to the target beam, and ends when the local contact indentation reaches to the maximum. The unloading phase begins at the local contact indentation reaching to the maximum, until the impactor totally departing from the geometrically imperfect GRC piezoelectric beam. According to the modified nonlinear Hertz contact proposed by Yang and Sun [50], the varying impact force F C t ð Þ during these two different phases can be predicted through two different equations.
During the loading phase [49], with the contact stiffness K C and the local contact indentation a t ð Þ where Y and m are the Young's modulus and Poisson ratio, respectively, and the subscript ''i'' and ''t'' respectively represent the impactor and the top surface of the geometrically imperfect GRC piezoelectric beam. Besides, w C denotes the deflection of the geometrically imperfect GRC piezoelectric beam at the impacted point.
During the unloading phase [49], where a max and F max are respectively the maximum local contact indentation and the corresponding impact force, and a 0 is the permanent indentation of the target beam. If the plastic deformation is not considered, a 0 ¼ 0.

Physical parameter model
It is assumed that the GRC beam is composed by a mixture of an isotropic matrix and GPLs with length l G , width w G and thickness h G . GPLs disperse into the GRC beam along its thickness direction. Both uniformly distributed (UD) and functionally graded (FG) patterns are taken into consideration, as shown in Fig. 1 f G denotes the total GPL weight fraction. q G and q M are respectively the mass density of the GPLs and matrix, and k ¼ 1; 2; Á Á Á ; N.
According to the Halpin-Tsai model [24] and the rule of mixture, the effective material properties, where the subscripts ''M'' and ''G'' represent matrix and GPLs, respectively. n LW and n WH are respectively the length-to-width ratio and width-to-thickness ratio of the GPLs, expressed by 3 Nonlinear dynamic equations In this study, the first-order shear deformation theory is used to estimate the deformation of the geometrically imperfect GRC piezoelectric beam. The axial and transverse displacementsũ x; z; t ð Þandw x; z; t ð Þof the novel beam can be presented by its mid-plane axial, transverse and rotary displacements u x; t ð Þ, w x; t ð Þ and u x; t ð Þ, combined with the initial transverse geometric imperfection Referring to von Kármán nonlinear displacementstrain relationship, the normal strain e xx and shear strain c xz of the geometrically imperfect GRC piezoelectric beam are gained as The linear constitutive relations for the kth GRC layer and piezoelectric layers can be expressed as where Y P and m P are Young's modulus and Poisson ratio of the piezoelectric layer, respectively, e 31 , k 33 and D Z are respectively the piezoelectric stress constant, permittivity constant and electric displacement of the piezoelectric layer. The electric potential variation is considered to be linear throughout the thickness of piezoelectric layer as Ref. [51], so the electric field E Z can be expressed as Based on Hamilton principle, the partial differential governing equations of motion for geometrically imperfect GRC piezoelectric beam under low-velocity impact are obtained as where the shear correction factor k s ¼ 5=6 and x C is the location of the impacted point. I j (j = 1, 2 and 3) are the inertia related terms The in-plane force N x , bending moment M x and shear forces Q x are respectively expressed by where the stiffness components A 11 , B 11 , D 11 and A 55 are respectively given by The thermal and electric related loads are written as Obviously, M P x ¼ 0. In this study, the geometrically imperfect GRC piezoelectric beam has hinged supported boundary conditions at both ends (x = 0, L). The corresponding boundary conditions are given as follows Introducing the following dimensionless quantities The partial differential governing Eqs. (23)- (25) can be rewritten in the dimensionless forms 4 Solution procedure In this paper, the differential quadrature (DQ) method [32,52] and Newmark-b method [53] combined with Newton-Raphson iterative process [13,54] are introduced to numerically analyze the nonlinear lowvelocity impact response of the geometrically imperfect GRC piezoelectric beam subjected to thermoelectro-mechanical loads. According to DQ method, the dimensionless displacement components U, W, w of the geometrically imperfect GRC piezoelectric beam and their rth derivatives with respect to X i are approximated in following forms where U m ; W m ; w m f g¼ U; W; w f g j X¼X m , l m X ð Þ is the Lagrange interpolation polynomial, C r ð Þ im is the weighting coefficients of the rth derivative at X ¼ X i , n is the total number of nodes distributed along the whole length of the GRC beam. In the present work, the distribution of nodes is given by the following equation Substituting Eqs. (39) and (40) into Eqs. (36)-(38), the dimensionless partial differential governing equations for nonlinear low-velocity impact response of geometrically imperfect GRC piezoelectric beam can be written as following ordinary differential equations With the DQ method, the boundary conditions in Eq. (34) can be rewritten as Hence, the above ordinary differential equations can be written as the following matrix equation M denotes the mass matrix, K L and K NL denote the structural linear and nonlinear stiffness matrix. Additionally, F and R are respectively time-dependent and time-independent load vectors, which respectively generated by impact and the geometrically imperfect combined with the thermo-electro-mechanical load Stiffness-proportional damping matrix C ¼ 2f=x l ð ÞK L [55][56][57] is considered to model energy dissipation of the impacted geometrically imperfect GRC piezoelectric beam, Eq. (44) can be written as where f is modal damping factor, and x l is the corresponding natural frequency. In this paper, x l is taken as the first natural frequency.
Assuming the initial conditions of the impactor as the governing equations of motion for the impactor, Eq. (2), and the ordinary differential equation of motion for the geometrically imperfect GRC piezoelectric beam, Eq. (49), can be solved parametrically by combining the Newmark-b method [53] with the Newton-Raphson iterative method [13,54] to analyze the time-dependent impact force F C t ð Þ and nonlinear dynamic behaviors of the low-velocity impact system. Figure 3 is the flowchart for the solution procedure, which is processing as below.
i. For the first time-step, assuming an initial impact force F C0 , the displacements y i1 and d 1 of the impactor and geometrically imperfect GRC piezoelectric beam can be respectively solved by Eq. (2) and Eq. (49) according to initial conditions d 0 , _ d 0 , € d 0 , y i0 and _ y i0 . ii. Applying the displacements obtained from step i into Eq. (3), a new impact force F C1 and local contact indentation can be predicted (Note that when the local contact indentation reaches the maximum, the new impact force needs to be calculated by Eq. (5)). If F C1 is close to F C0 , turn to step iii. If not, take the new impact force F C1 as the initial impact force, and repeat step i. iii. Treating the displacements and impact force obtained from steps i and ii as the initial conditions for the next time-step, repeat steps i and ii. iv. When the central displacement of the geometrically imperfect GRC piezoelectric beam remains unchanged, the iterative process ends.

Results and discussions
In this section, the nonlinear low-velocity impact behaviors are analyzed numerically for the GRC beam with geometric imperfections under thermo-electromechanical loads. The effects of weight fraction f G and distribution pattern of GPLs, mode and amplitude A 0 of imperfections, initial impact velocity V 0 , temperature variation DT and thickness ratio h P =h on the impact force and dynamic response of the system are discussed in detailed. Epoxy is selected as the matrix material of GRC beam, and the piezoelectric material is made by Polyvinylidene fluoride (PVDF). Material parameters [32,33,58] of epoxy, GPLs and PVDF are respectively m P ¼ 0:29; a P ¼ 145 Â 10 À6 /K; e 31 ¼ 0:05043 C/m 2 ; If there is no special instruction, length of the considered GRC beam L ¼ 50 mm, thicknesses of GRC host beam and piezoelectric layers are respectively h ¼ 10 mm and h p ¼ 0:5 mm, dimensionless amplitude of the geometric imperfection A 0 ¼ 0:02, the applied piezoelectric actuator voltage U 0 ¼ 400 V, temperature variation DT ¼ 50 K and damping factor f ¼ 0:1. The total number of the layers in the GRC host beam is N ¼ 12. The geometric parameters of GPLs are l G Â w G Â h G ¼ 2:5 lm Â 1:5 lm Â 1:5 nm and total weight fraction of GPLs is f G ¼ 0:3%. The impactor is made of steel, and its material parameters are Y i ¼ 207 GPa, q i ¼ 7960 kg/m 3 and m i ¼ 0:3. It is assumed that the impactor's radius is R i ¼ 5 mm, and the impactor with initial velocity V 0 ¼ 5 m/s impacts on the center (X C ¼ 0:5) of the top surface for the geometrically imperfect GRC piezoelectric beam.
Convergence and comparison studies are firstly performed to ensure the effectiveness and accuracy of the present method. Figure 4 calculates and compares the effect of geometric imperfection modes on (a) the highest impact force and the (b) maximum central deflections of the geometrically imperfect GRC piezoelectric beams with varying nodes number n during the impact process. Obviously, L-mode imperfections converge slowly. It is because L-mode imperfections have sudden changes in geometric characteristics (Fig. 2). Even G 1 -mode has the lowest impact force, it leads to the biggest center deflection. Hence, G 1 -mode imperfection is the most harmful geometric imperfection for the stability of the GRC piezoelectric beam under impact. In the following   analyses, we pay attention to the G 1 -mode imperfect GRC piezoelectric beam with nodes number n ¼ 13.
The comparison results are given in Table 1 and Fig. 5. In Table 1, the nonlinear frequency ratios are compared with Wu et al.'s results [59] for the intact and L 1 -mode imperfect FG-X CNTRC beams with nodes number n ¼ 39. Figure 5 compares the impact force history with Kiani et al.'s results [60] for the intact FG-X CNTRC beam subjected to a low-velocity impact with V 0 ¼ 3 m/s. The comparison results show that our results are in good agreement with the existing literature results. Figure 6 shows the influence of GPL distribution patterns (FG-X, UD, FG-O) on the impact force and impact response of the geometrically imperfect GRC piezoelectric beam. The case of pure epoxy piezoelectric beam's results is also given for comparing. GPLs distributing into the GRC core beam and the piezoelectric layers are respectively considered in Fig. 6a, b. It is known that GPL reinforcements can largely improve the stiffness of the matrix. Keeping in mind that the dynamic displacement and impact force F C are respectively closely related to the bending stiffness of the structure and contact stiffness of the impacted area. Therefore, when GPLs distribute into the GRC core beam, the varying GPL distribution pattern can significantly affect the dynamic behaviors, but has little effect on the impact force F C . To the Effect of GPL distribution patterns on the impact force and impact response of the geometrically imperfect GRC piezoelectric beam: a GPLs distributing into the GRC core beam and b GPLs distributing into the piezoelectric layers contrary, the impact force is distinctly affected by the different GPL distribution pattern when GPLs distribute into the piezoelectric layers. In addition, the FG-X pattern has the lowest central deflection W mid . For brief, only geometrically imperfect GRC piezoelectric beams with FG-X distributed GPLs are analyzed below. Figure 7 illustrates the influence of GPL weight fractions (f G ¼ 0:0%, 0.1%, 0.3% and 0.5%) on the (a) impact force history and (b) impact response of the geometrically imperfect FG-X GRC piezoelectric beams. f G ¼ 0:0% stands for the case of pure epoxy piezoelectric beam. As seen, both impact force F C and central deflection W mid decrease with the increase in GPL weight fractions, and the central deflection W mid is more sensitive to the varying weight fraction of GPLs. It indicates that even a small amount of GPLs can contribute to a great increase in bending stiffness of the geometrically imperfect FG-X GRC piezoelectric beam. Figure 8 displays the effect of the imperfection amplitude A 0 (A 0 = -0.02, -0.01, 0, 0.01, 0.02) on the (a) impact force and (b) impact response of the geometrically imperfect GRC piezoelectric beam. A 0 = 0 means the intact GRC piezoelectric beam. As seen, the impact force for the convex surface (A 0 \ 0) of the GRC beam impacted is obviously greater than that for the concave surface (A 0 [ 0) impacted. What's more, the impact force decreases with the increasing imperfection amplitude A 0 of the imperfection. Obviously, the final equilibrium position of (a) (b) Fig. 7 Effect of GPL weight fractions on the a impact force and b impact response for X-GRC piezoelectric beams (b) (a) Fig. 8 Effect of imperfection amplitude on the a impact force and b impact response of the geometrically imperfect X-GRC piezoelectric beam the geometrically imperfect GRC piezoelectric beam is also related to the sign of the imperfection amplitude A 0 . Bigger the absolute value of the imperfection amplitude is, further the final equilibrium position is apart from the z-axis. It is leaded by the direction and value of the component along z-axis of the force vector R in Eq. (49).
Note that the force vector R generated by the geometrically imperfect and the thermo-electro-mechanical load. R = 0 when the GRC piezoelectric beam is intact.
To illustrate the coupled effect of the geometric imperfection and thermo-electro-mechanical load on the nonlinear low-velocity impact behavior, Fig. 9 shows the effect of temperature variation (DT ¼ 0 K, 25 K, 50 K) on the nonlinear low-velocity impact response of (a) geometrically imperfect (A 0 ¼ À0:02 and 0.02) and (b) intact (A 0 ¼ 0) GRC piezoelectric beams. It can be observed that the impact force and impact response of the imperfect cases are much more sensitive to the varying temperature than those of the intact one. In other words, the influence of thermal load on the nonlinear low-velocity impact response is barely for the intact GRC piezoelectric beam, but can be largely amplified once the geometric imperfection appearing. Figure 10 demonstrates the (a) impact force history and (b) impact response of geometrically imperfect GRC piezoelectric beam impacted by different initial impactor velocities (V 0 ¼ 3 m/s, 5 m/s and 7 m/s). As seen that the increase in the initial impactor velocity V 0 leads to the remarkably increase in the maximum impact force, but obviously decreases the impacted time. Besides, the maximum center displacement increases a little for the increasing initial impactor velocity, which accompany with an increasing initial energy. Figure 11 demonstrates the effect of the thickness ratio (h P =h ¼ 0, 1/10, 1/5) of piezoelectric layer to GRC host beam on the (a) impact force F C and (b) central deflection W mid for the geometrically imperfect GRC piezoelectric beam. Here, we keep the thickness of GRC host beam unchanged. h P =h ¼ 0 represents GRC beam without any piezoelectric layers. Because the increase in piezoelectric layer thickness can improve the stiffness of the GRC piezoelectric beam, the impact force F C increases, and the absolute value of the central deflection W mid decreases.

Conclusions
The nonlinear low-velocity impact responses of GRC beam with geometric imperfection under thermoelectro-mechanical loads are investigated. The impact force is calculated by the modified nonlinear Hertz contact law. Based on the first-order shear deformation beam theory, von Kármán nonlinear displacementstrain relationship and Hamilton principle, the nonlinear dynamic equations are deduced, and solved by the DQ method and Newmark-b method combined (b) (a) Fig. 10 Effect of initial indenter velocity on the a impact force and b impact response of the geometrically imperfect X-GRC piezoelectric beam (b) (a) Fig. 11 Effect of thickness ratio on the a impact force and b impact response of geometrically imperfect X-GRC piezoelectric beam with Newton-Raphson iterative method. Comprehensive numerical results are presented to investigate the effects of the distribution pattern and weight fraction of GPLs, the mode and amplitude of imperfection, initial velocity of the spherical impactor, temperature variation, thickness of piezoelectric layer, as well as the coupled effect of geometric imperfection and thermo-electro-mechanical load on the low-velocity impact characteristics of the geometrically imperfect GRC piezoelectric beam. The major conclusions are given as follows.
(1) When calculating the impact behaviors of geometrically imperfect GRC piezoelectric beam via DQ method, the geometric imperfection with sudden changes require more nodes number. G 1 -mode is the most harmful imperfection for the stability of the GRC piezoelectric beam under impact. (2) When GPL reinforcements are distributed into GRC core beam, only FG-X pattern reduces the low-velocity impact response of geometrically imperfect GRC piezoelectric beam. For GPL reinforcements distributing into the piezoelectric layers, all the patterns are effective for reducing the central deflection of the geometrically imperfect GRC piezoelectric beam under impact. (3) The thermal load has barely effect on the impact behaviors of intact GRC piezoelectric beam. Once the geometric imperfection appearing, the thermal effect is largely amplified.