Nonlinear wave propagation in graphene incorporating second strain gradient theory

The exploration of graphene has attracted extensive interest owing to its significant structural and mechanical properties. In this research, we numerically investigate wave propagation in defect-free single-layer graphene, considering its geometrically nonlinear behavior through second-strain gradient elasticity. To capture the geometric nonlinearity, firstly, the nonlinear strain–displacement relations are introduced. The governing equation and associated boundary conditions are derived using Hamilton’s principle. Then, the weak form, including the element matrices, is established. The eigenvalue problem is solved for 2D wave propagation through periodic structures theory. Finally, the dynamical properties such as band structures, mode shapes, energy flow, and wave beaming effects are analyzed. The numerical results reveal that the geometric nonlinearity through the second strain gradient influences the wave propagation characteristics in graphene. The findings are significant and contribute to the understanding of graphene’s dynamic response, with implications for the engineering applications of graphene-based nanostructures


Introduction
As a two-dimensional nanostructure, graphene has gained significant attention due to its exceptional properties, including high strength, acoustic controlling capacity, and low energy consumption [1].Therefore, investigating the wave motion in graphene is of importance in gaining insight into its dynamic behavior, as it can be employed as a waveguide in various nano-and micro-electromechanical systems, including nanoelectronics and nanophotonics systems [2,3].Various experimental and theoretical investigations have been conducted to comprehend the mechanical and physical characteristics of graphene-based waveguides through diverse analysis techniques.
Among the initial approaches employed for analyzing graphene, atomistic methods such as the Molecular Dynamics (MD) simulation emerged prominently.For instance, Li et al. [4] reviewed the applications of MD simulation on mechanical properties of graphene sheet and introduced the capabilities of the method mentioned above on exploring mechanisms of graphene.Liu et al. [5] performed MD simulations to study the propagation of elastic waves in a nano-sized aluminum plate.Bo et al. [6] utilized the MD simulation to replicate the dispersion patterns of graphene within blends of polyvinylidene fluoride.
The results of their research demonstrate the accurate prediction of graphene distribution within the polystyrene phase.On the other hand, Arash et al. [7] developed a nonlocal elastic plate model to consider the scale effects of the wave propagation in graphene sheets.The result from finite element method is verified by MD simulations.Their findings indicate that the phase velocity decreases and approaches an asymptotic value as the width of the sheets reaches a sufficiently large size.Hu et al. [8] explored the transverse and torsional waves in both single-and double-walled carbon nanotubes (SWCNTs and DWCNTs), with a specific emphasis on assessing the impact of carbon nanotube microstructure on wave dispersion.Their findings revealed a good correlation between the wave dispersion predicted by nonlocal elastic cylindrical shell theory and the results obtained from MD simulations.Rahman et al. [9] proposed a MD approach to examine the thermal vibrations of graphene under various boundary conditions.The research was conducted on the examination of the natural frequency variations and their dependence on factors such as boundary conditions, length scale, and temperature.
However, the analysis of nanostructures such as graphene using atomistic methods is well-known to be computationally heavy and https://doi.org/10.1016/j.tws.2024.111713Received 6 December 2023; Received in revised form 13 February 2024; Accepted 13 February 2024 often impractical when dealing with large-scale structures comprising numerous molecules and atoms [10].In such cases, an alternative approach for analysis is to employ continuum mechanics.The mechanical behavior of nanostructures exhibits significant size effects and stiffness hardening phenomena [11].Nevertheless, the Classical Continuum Theory (CT), lacking the additional length scale parameters, fails to capture these special properties.Therefore, it is necessary to enhance the CT to adequately account for the size-dependent behavior, namely the non-local or long-range interactions [12], the large surfacevolume ratio [13], and the micro-deformations (e.g., micro-rotation and micro-stretch) [14], observed in nanostructures.
In order to interpret the properties of the nano-sized structures caused by the size effects, the non-classical continuum theories of elasticity have been proposed.Generally, these non-classical theories can be categorized into the First Strain Gradient (SG) theory [15], the Second Strain Gradient (SSG) theory [16], the couple stress theory [17], the modified couple stress theory [18,19], the micro-continuum theory [20,21], the surface elasticity theory [22], and the non-local elasticity theory [23].Notably, a singularity problem arises at the defect line of the double stresses by using the SG theory, in which the strain energy density comprises both strain and the first gradient of strain.To overcome this issue, the SSG theory was proposed, introducing the potential energy as a function of the first, second, and third gradient of displacement.The utilization of the SSG theory offers three significant advantages.Firstly, by considering double and triple stresses, it effectively eliminates the singularities present in physical fields, such as the elastic bend-twist tensor, which is singular in the SG theory [24], and the singularity of stress components at the crack tip [25].Secondly, it allows for the validation of higher-order deformations resulting from a higher-order characteristic length.Lastly, it demonstrates stiffness hardening due to non-local or long-range interactions.There seems to be no need to employ a third strain gradient theory since, according to the SSG theory, all physical state quantities are already smooth and non-singular.
Recently, researchers have conducted some explorations on the linear properties of nanostructures through the non-classical continuum theories.For example, Jalaei and Civalek [26] examined the dynamics of viscoelastic graphene sheet under periodic axial load via the SG theory.It is indicated that the effects of temperature reduce when the length scale parameter enhances.Ebrahimi et al. [27] reported a Navier-based analytical method and a semi analytical differential transform method to analyze the size-dependent vibration of nonlocal Euler-Bernoulli nano-beam.Liu et al. [28] studied the acoustic wave motion characteristics of the nanoplate through the modified couple stress theory.Zhang et al. [29] explored the mass exchanges with the bordering matrix and fracture, and reactive transport processes in the altered layer.What is more, Hu et al. [30] investigated the mechanical behavior of a Kirchhoff thin nanoplate based on the surface elasticity.It shows that there is a singularity near the rigid line inclusion tip for the bending moments.Zhang et al. [31] explained the dynamics of quadrilateral single-layered graphene sheets using the nonlocal elasticity theory.The result indicates that the skew angles and the magnetic field help increase the fundamental frequencies of single-layered graphene sheets.
Meanwhile, nanostructures, owing to their small size and unique properties, exhibit nonlinear behavior [32].This behavior can influence the performance of structures in wave propagation, including frequency shifting and consequently, a shift in the location of the band gap [33].Understanding nonlinear wave propagation is crucial for designing and optimizing the performance of nanostructures.Extensive research has been conducted on the nonlinear characteristics of macro-sized structures [34][35][36][37].However, there is a scarcity of literature concerning the nonlinear dynamical behaviors of nanostructures.Karimiasl et al. [38] investigated nonlinear vibration behavior of multiscale nanocomposites nanoshell.Esfahani et al. [39] studied the nonlinear behavior of a nonlocal strain gradient functionally graded nano resonator.It is shown that the surface effects induce the stiffness hardening behavior.Behrouz et al. [40] proposed a new biological layer detection strategy based on the nonlinear couple stress theory and the results were validated qualitatively with previous literature.Sahmani et al. [41] captured the size dependencies in the nonlinear nano-beams by the nonlocal SG theory.It shows that the increment from the size dependency is more than the reduction caused by the nonlocality.On the other hand, Audoly and Lestringant [42] proposed an approach to calculate the equivalent one-dimensional nonlinear models of slender structures.Khayat et al. [43] analyzed the nonlinear dynamic properties of a sandwich structure reinforced with the graphene.According to the results, the sensitivities of the responses have the highest values related to the uncertainty sources.Furthermore, various experimental explorations on the mechanical properties of graphene have been conducted.For instance, Xu et al. [44] experimentally explored the interfacial mechanical characteristics of a sizable monolayer graphene adhered to a flexible Polyethylene Terephthalate (PET) substrate.It shows that the mechanical characteristics of the interaction between a sizable graphene sheet and a substrate are influenced by various factors, including the dimensions of the graphene, the texture of the substrate material, and the roughness of the substrate.Hadden et al. [45] studied the mechanical properties of Graphene Nanoplatelets (GNPs) through multiscale modeling and experiments.The result indicates that the accuracy of the multiscale modeling approach offers valuable physical insights into the mechanical behavior of the composite.On the other hand, Lopez et al. [46] summarized three experiments for the investigation on the effect of rippling on the mechanical properties of graphene.The consistent findings from the experiments suggest that the mechanical properties of graphene are affected solely by rippling with wavelengths of some certain ranges.Tu et al. [47] conducted experimental analysis on the interfacial mechanical properties of graphene on Self-Assembled Monolayers (SAM).Their findings indicate that the head group chemistry of a SAM can notably influence the out-of-plane elastic modulus of the graphene-SAM heterostructure.
On the other hand, the periodic nanostructures with special properties have been widely applied in various engineering fields.In order to investigate the dynamics of these structures mentioned above, a Wave Finite Element Method (WFEM) [48] has been put forward.This approach represents an advancement from the conventional Finite Element (FE) method and allows for modeling of complex structures.By applying the principle of periodic structures theory, such as the Floquet-Bloch theory [49], the WFEM can simplify a periodic structure into a single unit cell, reducing the computational burden associated with the analysis.What is more, the homogenization technique [50], including the asymptotic homogenization method [51] and equivalent strain energy method [52], has also been widely used to investigate the periodic structures.
The authors have previously introduced adaptations of the WFEM framework for the analysis of linear 1D bending and torsion wave propagation [53], as well as linear 2D wave motion in periodic structures [54], within the context of strain gradient theories.In this study, for the first time to the author's knowledge, WFEM is used to investigate the nonlinear 2D wave propagation in graphene within the framework of SSG theory.The integration of the SSG theory with WFEM holds importance as it enables the explanation of size-dependent characteristics of micro-sized structures using the SSG theory and facilitates the exploration of the dynamical properties of complex periodic structures through WFEM.Furthermore, to capture the geometric nonlinearity, the nonlinear strong and weak forms of graphene are derived based on the SSG theory for the first time.Lastly, this paper gives new insights regarding the properties observed in graphene, which have considerable implications for the outcomes but cannot be explained by the classical theories.The aforementioned novel findings in this paper enhance the understanding of nonlinear dynamic properties of graphene and make contributions to the practical implementation of graphene-based nanostructures in engineering applications.The main objectives of this work are, firstly, to introduce the nonlinear strong and weak forms including the element matrices through the SSG theory, and, secondly, to investigate the nonlinear 2D wave motion properties, such as band structures, mode shapes, energy flow, and wave beaming effects, by solving the eigenvalue problem through the WFEM.
The paper's structure is the following: in Section 2, the geometrically nonlinear mass-spring model and element matrices of graphene are confirmed through the SSG theory.Afterwards, in Section 3, the free nonlinear 2D wave propagation formulations are calculated by solving the eigenvalue problem within the WFEM framework.In Section 4, the dynamical properties of graphene are explored.Finally, some conclusions are presented in Section 5.

Graphene modeling
In this part, firstly, the geometrically nonlinear mass spring model of graphene is introduced.Then, the governing equation and associated boundary conditions are confirmed by the SSG theory.Finally, the weak form including the element matrices is calculated.

Nonlinear mass spring model
As shown in Fig. 1(a), a nano-electromechanical resonator consisting of a graphene-based waveguide and electrical signal ports is proposed in this paper.The intensity at electrical excitation port is modulated at frequency , causing a periodic motion.The motion can be detected by monitoring the reflected light intensity at electrical detection port.As shown in Fig. 1(b), the carbon atoms interact through out-of-plane -bonds and in-plane -bonds.In this paper, the graphene is assumed to be a defect-free single layer, referring to a structurally perfect form of graphene without any defects or imperfections in its carbon lattice.The in-plane -bonds are considered only.Here, we focus on the continuum mechanics of the graphene.As addressed in Fig. 1(c), a geometrically nonlinear mass spring model is used to reflect the interaction between two in-plane carbon atoms.

Finite element discretization
In order to confirm the mechanical model of the graphene, firstly, the strain energy density  via the SSG theory [16] is introduced, as below: where  = (∇) indicates the first gradient of displacement,  = ∇∇ represents the second gradient of displacement, and  = ∇∇∇ denotes the third gradient of displacement. is the classical stiffness tensor determined by the conventional  é parameters  and .7) , and  = (  ) (=1,…,3) are the nonclassical stiffness tensors associated with higher-order parameters   ,   , and   , respectively. 0 =  0 ( 0 ) is the dimension of force related to a cohesion modulus  0 .The parameters   , capturing variations in strain throughout the material, are related to the first-order spatial derivatives of the strain tensor.The parameters   , elucidating the curvature of the strain field and offering insights into the spatial variations of strain, describe the second-order spatial derivatives of the strain tensor.Furthermore, the parameters   encompass the coupled influence of the strain tensor and the second-order spatial derivatives of the strain tensor.The higher-order parameters mentioned above can be confirmed by the inter-atom potential function method [55].Here, ∇ signifies the gradient operator, and  denotes the displacement vector.In this study, as shown in Fig. 1(c), the graphene is assumed as a defectfree single layer.Consequently, only the in-plane -bonds exhibiting tension behavior along the  direction are taken into account, while no displacements are considered along the  and  directions.The displacement vector can be expressed as: in which (, ), representing the displacement along the local coordinate , is considered only.The Green-Lagrange strain [56] and its first and second gradient, reflecting the geometric nonlinearity, are of the form: ) 2 . ( Substituting Eq. (3) into Eq.( 1), as a result, the re-expression of the strain energy density within the SSG theory is given by the following equation: where )∕ = ∕, and   = ( 1 + 2 + 3 )∕ = ∕ are the higher-order material parameters existing in the SSG theory. and  are related to the Young's modulus , the Poisson's ratio , and the shear modulus , as Here,  =  for 1D mass spring tension model.Then, the constitutive equations can be written as: Integrating the strain energy density in Eq. ( 4) over its volume to obtain the strain potential energy, one arrives: in which  is the cross-sectional area and   denotes the length of the structure along .The definitions of  ,   ,   ,   , and  0 are addressed in Appendix A. On the other hand, the kinetic energy is represented by the combination of the classical component and the non-classical component, and it can be formulated in the form: where  1 and  2 are the higher-order inertia parameters.Furthermore, the work done by external classical force and higher-order forces are given by: in which  (, ) means the distributed force,  0 denotes the classical nodal force, and  1,2 are the higher-order forces, which are assumed to be applied only at the both end sections.Next, in order to illustrate the strong form, the Hamilton principle is used via the SSG theory as follows: Where  ,  , and  represent the variations in kinetic energy, strain energy, and external work done, respectively.Then, by substituting Eqs. ( 6), (7), and (8) into Eq.( 9) and performing mathematical calculations using the variation method, the governing equation can take the following form: In order to guarantee the continuity of high-order displacement, the six-term polynomial function is used to interpolate the scalar field (, ) inside the 1D element: Substituting Eq. ( 12) into Eq.( 11), the nodal displacement vector   () can be re-written as: where  1 = {1, 0, 0, 0, 0, 0},  2 = {0, 1, 0, 0, 0, 0},  3 = {0, 0, 2, 0, 0, 0},  4 = {1,   ,   12) with Eq. ( 13), the displacement vector can be illustrated by using the  2 continuum Hermite interpolation function and nodal displacement vector: in which the interpolating function () takes the following form: where N 0 1 (), N 1 1 (), N 2 1 (), N 0 2 (), N 1 2 (), and N 2 2 () can be written as: Multiplying   (, ) as a test function for the governing equation, and performing the part-by-part integration respect to .Then, by utilizing Eq. ( 14) and   (, ) = ()   (), the geometrically nonlinear governing equation of an element in matrix form can be confirmed as: where in which   and   represent the linear stiffness matrix and mass matrix, respectively.  1 and   2 indicate the nonlinear stiffness matrices.  denotes the linear force vector.Superscript ( ′ ) is the partial derivative with respect to coordinate .

Nonlinear wave motion
After obtaining the nonlinear element matrices and force vector, in this part, firstly, the dynamic equilibrium equation of a unit cell is deduced.The frequency and displacement components for the perturbed system are defined.Then, the 2D wave motion properties are analyzed by solving the eigenvalue problem through the periodic structures theory.

Nonlinear dynamic formulation of a unit cell
According to the WFEM, only one unit cell or substructure is modeled through the traditional finite elements.The dynamic equilibrium formulation of a unit cell can be written as: where ,  1 ,  2 , and  are the unit cell stiffness and mass matrices assembled by   ,   1 ,   2 , and   . = ∕ is defined as a damping matrix, in which  is the damping loss factor caused by internal friction within the material. is nodal displacement vector. represents nodal force vector.A scheme for the 2D periodic structure and its unit cell is given in Fig. 2(a).The nodal DOFs in the unit cell are divided into six in which  is the frequency and  indicates the time.To characterize the perturbed response of the structural system, a perturbation parameter  is introduced.Using the Linstedt-Poincaré method [57], Eq. ( 19) can be asymptotically expanded.The nonlinear stiffness matrices can be written as: where perturbation parameter  is used to quantify the magnitude of the nonlinearity.The frequency and displacement components of the perturbed system can be determined using the following first-order Linstedt-Poincar é expansion: Substituting Eq. ( 20), Eq. ( 21), and Eq. ( 22) into Eq.( 19), the linear contribution •( 0 ) and the first-order contribution •( 1 ) can be separated as: Furthermore, the harmonic displacement and force components can be expressed in the form: where û0 = Q 0 Φ 0 and F = Q F Φ F , with the amplitude  0∕F and the eigenvector Φ 0∕F .Then, for the •( 0 ), the dynamic stiffness matrix can be re-expressed as  = K −  2 0  in frequency ( 0 ) domain with K = (1 + i), as a result: It is important to highlight that the internal DOFs remain unaffected by external loads, as the coupling actions are confined to the system's boundaries exclusively.So,  I = 0 and Eq. ( 25) will be rewritten as: where D 11 =  Bd,Bd , D 12 =  Bd,I , D 21 =  I,Bd , and D 22 =  I,I .The expression of DOFs in accordance with the periodic structures theory involves the utilization of propagation constants   and   [48]: in which   =  −i    and   =  −i    .The values of wavenumbers   and   vary within the boundaries of the first Brillouin zone, namely [−∕  , ∕  ] and [−∕  , ∕  ], respectively, as depicted in Fig. 2(b).Meanwhile, in the case of free wave propagation, when the sum of nodal forces exerted on all elements connected to nodes c 1 , L, and B are zero, one arrives:

Corrected frequency solution
In this part, the 2D wave propagation problem will be discussed.The nodal DOFs of a unit cell can be expressed as: with where ,  s ,  sn ,  sm , and  i indicate the identity matrix of size s, sn, sm, and i, respectively.On the other hand, the nodal forces can take the form: with Then, combining Eq. ( 26), Eq. ( 29), and Eq. ( 31), this yields: To solve Eq. ( 33), we follow the approach of fixing the values of   and   , enabling the calculation of the corresponding  0 values.These  0 values are then sorted in ascending order as  , 0(1,2,3,…,,…,) , where  represents the th term of wavenumber along the -direction,  denotes the th term of wavenumber along the -direction,  is the th term of frequency under  and .The formation of the th slowness surface is established based on  , 0() .At a given point (  ,   ), the energy flow vector, also known as the Poynting vector, is equal to the gradient of the slowness surface.When considering the iso-frequency contour of the slowness surface, the Poynting vector is perpendicular to the contour curves.A convenient way to represent the wave propagation is through the band structure, which can be obtained by plotting the wavenumbers along the contour O-A-B-C-O, as illustrated in Fig. 2(b).
On the other hand, submitting Eq. ( 29) into the part of  0 in Eq. ( 24), and doing the same calculation on  1 , one arrives: Then, submitting Eq. ( 34) into the part of •( 1 ) in Eq. ( 23), as a result: Here, it should be pointed out that the expression  i was equivalently replaced by cos, cos 3  was linearized as 3  4 cos.Due to its disruptive effect on the symmetry property, the stiffness matrix  1 was disregarded [58].Meanwhile, the eigenvector of the left side of Eq. ( 35) is identical to the one of Eq. ( 33), namely ûb So, the Eq. ( 35) will be re-expressed as: The corrected frequency  1 can be calculated by submitting  2 in Eq. ( 21) into Eq.( 36).As a result, frequency  can be subsequently determined by combining the solutions of the linear zero-order system and the perturbed system.This results in the following expression: On the other hand, according to Eqs. ( 22) and ( 34), the nonlinear eigenvector can be illustrated by determining  = QΦ i , as follows: and the nonlinear forced response can be expressed as: where  is the amplitude of displacement vector of global structure. * = K * −  2  * with K * = (1 + i) * .The global stiffness  * and mass  * matrices are assembled by the unit stiffness () and mass () matrices, respectively. * is the amplitude of global force vector.

Numerical results
The numerical exploration of graphene's nonlinear characteristics using the SSG theory holds significant value, as it provides essential insights for the practical utilization of graphene-based nanostructures in engineering applications.In this part, the dynamic behaviors such as band structure, mode shapes, energy flow vector fields, and wave beaming effects are illustrated.

2D wave propagation
In this section, the band structures of graphene are introduced to investigate the 2D wave propagation.Fig. 4(a) illustrates the first ten branches of the normalized frequency spectrum along the boundary O-A-B-C-O of the irreducible first Brillouin zone.The results obtained from the linear SSG theory are represented by the continuum lines, while the dotted lines correspond to the linear CT results.It shows that at low frequencies, the linear SSG theory closely aligns with the linear CT curve.However, as the frequency increases, the disparities between the linear SSG theory and linear CT results become more pronounced.Furthermore, the frequency values obtained from the linear SSG theory are higher than those from the linear CT at the same -space position.This discrepancy can be attributed to the potential energy density in the linear SSG theory with the additional first gradient of strain and second gradient of strain.The dynamical equilibrium equation by linear SSG theory becomes a high-order partial differential function encompassing both classical and non-classical components.The inclusion of non-classical components involving size effects causes the eigenvalue calculated from the linear SSG theory to surpass the value obtained from the linear CT at the same -space position, and leads the stiffness to be increased.To validate the band structure derived from the linear SSG, a linear Component Mode Synthesis (CMS) method [63] is employed.As depicted in Fig. 4(a), the CMS results match well with the WFEM results, especially at low frequencies.Additionally, when the higher-order parameters are set to zero, the linear SSG theory converges to the linear CT, providing a reasonable way to verify the accuracy of the simulation results obtained using the linear SSG theory.
On the other hand, Fig. 4(b) exhibits the band structures of the first five branches obtained through both linear and nonlinear SSG theories.A comparison between the two reveals that the results from linear SSG theory are close to the ones from nonlinear SSG theory for the first and fifth branches.However, along the O-A-B-C-O contour of the first Brillouin zone, the frequency values obtained from the nonlinear SSG theory surpass those from the linear SSG theory for the second branch.Additionally, concerning the third branch, the frequency values obtained from the nonlinear SSG theory are higher along the O-A and B-C-O contours but nearly equivalent to those obtained from the linear SSG theory along the A-B contour.Regarding the fourth branch, the frequency values from the nonlinear SSG theory are higher along the O-A-B contour but approach those from the linear SSG theory along the B-C-O contour.The difference between linear and nonlinear outcomes can be elucidated by considering the potential energy density in the nonlinear SSG theory.This density is a function of strain, the first gradient of strain, and the second gradient of strain, which includes a geometric nonlinear component.Consequently, the dynamical equilibrium equation becomes a high-order partial differential function comprising both linear and nonlinear parts.The incorporation of a nonlinear element in the strain energy density leads to higher frequency values at the same position in -space, resulting in stiffness hardening of the material.This phenomenon illustrates the influence of   nonlinearity on structural properties.To validate the nonlinear results, the nonlinear CMS is employed, and its findings match well with the nonlinear WFEM results.

Impact of higher-order parameters
As discussed in the previous section, the band structures within the SSG theory framework are affected by the higher-order parameters   ,   ,   ,  1 , and  2 .To examine the influence of these parameters on the outcomes, each parameter is multiplied by a ratio , and the effects are analyzed by varying the value of  from 0 to 5 in increments of 0.5.There are 5 different cases:   changes only for case 1,   changes only for case 2,   changes only for case 3,  1 changes only for case 4, and  2 changes only for case 5.In this study, we focus on analyzing the first normalized frequency at the -space position A 1 .Fig. 5 presents the impact of the higher-order parameters on the results.Here, the frequency growth rate   is defined as lg|( (+1) −  () )|, where  as number of  ranges from 1 to 10. Fig. 5(a) addresses the influence of the higher-order material parameters (  ,   , and   ).The figure demonstrates that as the values of these parameters increase, the frequency growth rate also increases, except for the nonlinear   , which initially decreases and then increases.For   and   , the frequency growth rate considering nonlinearity is smaller than the linear one, except for  = 5 in the case of   .For   , the frequency growth rate considering nonlinearity is greater than the linear one at high  value.Furthermore, Fig. 5(b) illustrates the impact of the higher-order inertia parameters ( 1 and  2 ).The results indicate that as the value of  increase, the frequency growth rate initially increases and then decreases.At the same  value, the frequency growth rate considering nonlinearity is higher than the linear rate.Notably, the higher-order material parameter   exhibits the most significant influence on the wave propagation characteristics.
To specifically discuss the impact of higher-order parameters on the results, the normalized frequencies at A 1 using both linear and nonlinear SSG theories for various  values are calculated.As presented in Table 1, within the SSG theory framework, the nonlinear case consistently produces higher frequency values than the linear case, indicating a significant influence of nonlinearity on wave propagation.Furthermore, for case 1 (effect of   ), the frequency value obtained from both the linear and nonlinear SSG theory increases as  increases.For case 2 (effect of   ), the frequency value obtained from both the linear and nonlinear SSG theory decreases firstly and then increases as  increases.For case 3-5, the influences of   ,  1 , and  2 continue to diminish as  increases.Notably, the impact of   is more pronounced than that of other higher-order parameters when  exceeds 1.5.

Mode shapes of a unit cell
The mode shapes provide insight into the propagation of waves in graphene.The normalized mode shapes at the intersecting points,  ).Upon normalization, the nonlinear eigenvector will be same as the linear eigenvector.For this discussion, only the normalized mode shapes are considered.Figs.6(a) -6(d) present the mode shapes corresponding to wavenumbers   = ∕  and   = 0. Figs.6(e) -6(h) display the mode shapes when the wavenumbers are   = 0 and   = ∕  .These figures reveal that the mode shapes of the hexagonal carbon ring at -space position A 3 exhibit similarities to those at position C 4 .However, the mode shapes of the left and right branches connected to the hexagonal carbon ring at position A 3 differ from those at position C 4 .Additionally, at the same -space position, the mode shapes undergo changes as the frequency increases.For example, at -space positions A 1 and A 2 , the mode shapes at the junctions between the hexagonal carbon ring and the two branches appear smooth.As the frequency increases, such as at space positions A 3 and A 4 , the mode shapes at these junctions become non-smooth.Moreover, the simulation results indicate that the wave has a greater influence on the hexagonal carbon ring than on the two branches at both low and high frequencies.

Energy flow vector fields
flow fields by linear and nonlinear SSG theories related to the first (low) and fifth (high) frequencies in the first Brillouin zone are studied in this part.Fig. 7 illustrates the energy flow, where the arrow's direction indicates the energy flow direction and the arrow's length represents the gradient value of the energy flow.At low frequencies (see Figs. 7(a) and 7(b)), the behavior of graphene resembles that of a homogeneous plate.The dynamic energy spreads in all directions and is perpendicular to the iso-frequency contour.From the center point to the boundary, the change of energy gradually decreases as the frequency rises, becoming more concentrated.The structure maintains its stability in all directions.However, at higher frequencies, such as 31.6 (Figs.7(c) and 7(d)), the majority of the energy is limited to the  direction, resulting in stop bands in the  direction.The high-energy waves, which cause significant vibrations in the structure, predominantly propagate in the  direction.Furthermore, at low frequencies, the iso-frequency contour and wave propagation behavior of the linear model closely resemble those of the nonlinear model.However, at high frequencies, a noticeable discrepancy emerges between the linear and nonlinear models.The wave propagation range in the nonlinear model is wider than that in the linear model.
The propagation of plane waves can be depicted by measuring the distance between any point on the iso-frequency contour and the central point of the figure.At lower frequencies, the iso-frequency contour by nonlinear SSG theory lies within that by linear SSG theory, signifying a smaller wavenumber from nonlinear SSG theory compared to linear SSG theory at the same wave propagation angle.In contrast, at higher frequencies, the iso-frequency contour by nonlinear SSG theory extends beyond that by linear SSG theory, indicating a larger wavenumber from nonlinear SSG theory compared to linear SSG theory.

Wave beaming effects
In this part, the wave beaming effects under Born-von Karman boundary conditions [63] are investigated by exploring the harmonic displacement fields.In damped systems, the response of harmonic displacement is treated independently of the boundary conditions, as the damping loss factor  limits the input power from external forces, resulting in minimal energy reaching the system boundary.Consequently, the system's response closely approximates that of an infinite system.In this work, the graphene is constructed by 20 unit cells along  direction and 20 unit cells along  direction with free boundary condition.The   loss factor  = 0.2 is used and a harmonic force with unit amplitude at the central location of the graphene along  direction is applied.
The displacement amplitudes, normalized with respect to the amplitudes at the central point, are depicted in Fig. 8.The linear and nonlinear SSG theories are employed to generate displacement fields.The directional patterns of the harmonic displacement response match the predictions derived from the energy flow vector fields.Moreover, the linear SSG theory accurately predicts the wave propagation range at low frequencies, matching the results obtained from the nonlinear SSG.However, discrepancies arise at high frequencies, where the wave propagation range estimated by the nonlinear SSG theory differs from that predicted by the linear SSG theory.Notably, the wave propagation range estimated by the nonlinear SSG theory surpasses that of the linear SSG theory when the normalized frequency reaches 31.6.
Lastly, Fig. 9 depicts the iso-displacement contours calculated using both linear and nonlinear SSG theories.These contours are obtained under a normalized displacement w = 0.2.The normalized frequency was set to 31.6, representing a high frequency.Both linear and nonlinear SSG theories demonstrate wave spreading in the  direction, indicating strong and evident beaming in that direction.The presence of geometric nonlinearity in the structure gives rise to distinct properties, leading to a change in the range of wave propagation.When employing nonlinear SSG theory at the same frequency, the wave propagation range surpasses that by linear SSG theory.Furthermore, as one moves from the midpoint of the graphene to the edge along the  direction, the wave propagation range initially diminishes and then expands by nonlinear SSG theory.In contrast, the wave propagation range shows a consistent increase by linear SSG theory.

Conclusions
This paper presents an analysis of 2D wave propagation in geometrically nonlinear graphene waveguides.Based on the results of this study, some conclusions are presented and listed as follows: (1) The study takes into account the higher-order strain gradient effect.To carry the numerical modeling, the strain-displacement relations considering geometric nonlinearity based on the SSG theory are introduced.The governing equation and boundary conditions are derived using Hamilton's principle.The weak form, which includes the element matrices, is established.
(2) Subsequently, the band structures of graphene are analyzed.The investigation compares the frequency values obtained from the SSG theory with those from the CT at corresponding -space positions.Remarkably, the SSG theory yields higher frequency values than the CT at the same -space position and leads to stiffness hardening phenomena in nanostructures.The different outcomes between the linear and nonlinear results arise due to the incorporation of geometric nonlinearity.The consideration of a nonlinear component in the strain energy density contributes to elevated frequency values in the nonlinear results when compared to the linear results.The influence of higher-order parameters on the band structures is demonstrated.Notably, the higherorder material parameter   exhibits the most significant impact on the results.
(3) Furthermore, an analysis is conducted on the mode shapes, which offer insights into wave propagation in graphene.The findings indicate that the mode shapes of the hexagonal carbon ring located at -space position A exhibit similarities to those at -space position C.However, the mode shapes of the left and right branches connected to the hexagonal carbon ring at -space position A differ from those at -space position C.The study also investigates the energy flow vector fields using linear and nonlinear SSG theories associated with the low and high frequencies.At low frequencies, graphene behaves similarly to a homogeneous plate.The iso-frequency contour and wave propagation behavior of the linear model closely resemble those of the nonlinear model.At high frequencies, a significant portion of the energy is confined to a single direction.The wave propagation range in the nonlinear model is wider compared to that in the linear model.
(4) Lastly, the wave beaming effects are investigated by exploring the harmonic displacement fields through the linear and nonlinear SSG theories.The observed directional patterns of the harmonic displacement response align with the predictions derived from the energy flow vector fields.At low frequencies, the wave propagation range by the linear SSG theory is consistent with the one by the nonlinear SSG theory.However, discrepancies arise at high frequencies, where the estimated wave propagation range from the nonlinear SSG theory deviates from the prediction made by the linear SSG theory.To further illustrate this, the iso-displacement contours were calculated.Both the linear and nonlinear SSG models exhibit wave spreading in the  direction, providing strong beaming effects.The findings presented in this study can enhance the comprehension of the dynamic response of graphene, with potential implications for engineering applications involving graphene-based nanostructures.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.The design and structure of a nano-electromechanical resonator.(a): The schematic of a graphene resonator which is constructed by electrical excitation port, electrical detection port, electrical tuning, and on-chip periodic graphene-based waveguide.(b): Two interacting carbon atoms in graphene.(c): The nonlinear mass spring model.

Fig. 2 .
Fig. 2. A unit cell in 2D periodic structures and the first Brillouin zone.L, R, B, T denote left, right, bottom, top boundaries DOFs.I means internal DOFs. 1 ,  2 ,  3 and  4 indicate the DOFs on four corners.  and   are the length of the unit cell along the  and  directions, respectively.

Fig. 3 .
Fig. 3.A defect-free single layer armchair graphene sheet and its unit cells.(a): The graphene is constructed by   unit cells along  and   unit cells along , respectively.(b): A unit cell and its geometric description.

B
.Yang et al.

Fig. 4 .
Fig. 4. The band structures of the defect-free single layer armchair graphene.(a): The first ten branches of band structures through linear SSG theory and CT.(b): The first five branches of band structures via linear and nonlinear SSG theory.  =  0 indicates linear,   =  means nonlinear.

Fig. 5 .
Fig.5.The impact of higher-order material and inertia parameters on the first frequency at -space position A 1 of band structures.(a): The influence of higher-order material parameters on the frequency growth rate.(b): The influence of higher-order inertia parameters on the frequency growth rate.

Fig. 6 .
Fig. 6.The first four mode shapes of a unit cell by SSG.(a): The mode shape at A 1 .(b): The mode shape at A 2 .(c): The mode shape at A 3 .(d): The mode shape at A 4 .(e): The mode shape at C 1 .(f): The mode shape at C 2 .(g): The mode shape at C 3 .(h): The mode shape at C 4 .namely A 1 , A 2 , A 3 , A 4 , C 1 , C 2 , C 3 , and C 4 as shown in Fig. 4(b), are depicted in Fig. 6.It should be noted that in Eq. (37), the nonlinear eigenvector is constructed by the nonlinear amplitude ((1 + )Q b 0 ∕Q) and the linear eigenvector (Λ R Φ b 0 ).Upon normalization, the nonlinear eigenvector will be same as the linear eigenvector.For this discussion, only the normalized mode shapes are considered.Figs.6(a) -6(d) present the mode shapes corresponding to wavenumbers   = ∕  and   = 0. Figs.6(e) -6(h) display the mode shapes when the wavenumbers are   = 0 and   = ∕  .These figures reveal that the mode shapes of the hexagonal carbon ring at -space position A 3 exhibit similarities to those at position C 4 .However, the mode shapes of the left and right branches connected to the hexagonal carbon ring at position A 3 differ from those at position C 4 .Additionally, at the same -space position, the mode shapes undergo changes as the frequency increases.For example, at -space positions A 1 and A 2 , the mode shapes at the junctions between the hexagonal carbon ring and the two branches appear smooth.As the frequency increases, such as at space positions A 3 and A 4 , the mode shapes at these junctions become non-smooth.Moreover, the simulation results indicate that the wave has a greater influence on the hexagonal carbon ring than on the two branches at both low and high frequencies.

B
.Yang et al.

Fig. 7 .
Fig. 7. Energy flow vector fields by SSG.(a): Linear energy flow related to the first frequency in the first Brillouin zone.(b): Nonlinear energy flow related to the first frequency in the first Brillouin zone.(c): Linear energy flow related to the fifth frequency in the first Brillouin zone.(d): Nonlinear energy flow related to the fifth frequency in the first Brillouin zone.The yellow line indicates the iso-frequency contour.

Fig. 9 .
Fig. 9. Normalized forced response contours by SSG theory under w=0.2.  denotes the number of unit cells along .  means the number of unit cells along .
.The total DOFs in a 1D element can be written as: with associated boundary conditions presented in Appendix B. The detailed expressions of  1 ,  2 ,  3 , and  4 are addressed in Appendix C. The subsequent stage involves converting the strong form into its corresponding weak form.Based on the SSG theory, there are three DOFs,   ,   2 ,  = 1, 2, on each node of an element with length s    s      s

Table 1
The normalized linear and nonlinear frequencies at A 1 via different cases.