Namnabat et al. (cf., [Carbon Letters, https://doi.org/10.1007/s42823-020-00194-2]) employ the classical approach of Li and Chou (cf., [Int J Solids Struct 40: 2487–2499]) to the implementation of the molecular structural mechanics method using the Bernoulli–Euler beam elements for nonlinear buckling analysis of double-layered graphene nanoribbons. However, more recent studies by Eberhardt and Wallmersperger (cf., [Carbon 95: 166–180]) and others (see, e.g., [Int J Eng Sci 133: 109–131]) have shown that the classical approach of Li and Chou poorly reproduces both in-plane and out-of-plane mechanical moduli of graphene. We have shown that the 2D beam-based hexagonal material used by Namnabat et al. poorly simulates the mechanical moduli of graphene, especially the bending rigidity modulus, and this material cannot be used for the buckling simulation of graphene sheets (or nanoribbons). In addition, it is noted that in Int J Eng Sci 133: 109–131, a modification of the classical approach of Li and Chou is given which exactly reproduces both in-plane (2D Young’s modulus and Poisson’s ratio) and out-of-plane (bending rigidity modulus) mechanical moduli of graphene using beam elements.

Double-layer graphene nanoribbons promise potential application in nanoelectromechanical systems and optoelectronic devices, and knowledge about mechanical stability is a crucial parameter to flourish the application of these materials at the next generation of nanodevices. In this paper, molecular mechanics is utilized to investigate nonlinear buckling behavior, critical buckling stress, and lateral deflection of double-layered graphene nanoribbons under various configurations of stacking mode and chirality. The implicit arc-length iterative method (modified Riks method) with Ramm’s algorithm is utilized to analyze the nonlinear structural stability problem. The covalent bonds are modeled using three-dimensional beam elements in which elastic moduli are calculated based on molecular structural mechanics technique, and the interlayer van der Waals (vdW) interactions are modeled with nonlinear truss elements. An analytical expression for Young’s modulus of nonlinear truss elements is derived based on the Lennard–Jones potential function and implemented in numerical simulation with a UMAT subroutine based on FORTRAN code to capture the nonlinearity of the vdW interactions during the buckling analysis. The results indicate that the highest critical buckling stress and the minimum lateral deflection occur for armchair and zigzag chirality, both with AB stacking mode, respectively. Moreover, the critical buckling stress is found to be directly dependent on the mode shape number regardless of in-phase or anti-phase deflection direction of layers. Lateral deflection exhibits a similar trend with mode shape in anti-phase mode; however, it is decreasing by increasing mode shape number in in-phase mode.