Abstract: Within the C⁰-higher-order shear deformation theory (C⁰-HSDT) framework with seven variables, the classical third-order shear deformation theory (TSDT), sinusoidal shear deformation theory (SSDT) and exponential shear deformation theory (ESDT) are selected. Based on the von Karman theory of large deflection, a finite element model is established to analyze the large deflection bending behavior of functionally graded graphene-reinforced composite (FG-GRC) plates. The equivalent Young's modulus and equivalent Poisson's ratio of the FG-GRC plate are determined using the modified Halpin-Tsai micromechanical model and the mixing rule. The model satisfies the zero traction condition on the top and bottom surfaces of the plate, which avoids the difficulties in constructing C¹-continuous elements in the finite element method (FEM). To overcome errors introduced by the artificial rotational variables in traditional C⁰-HSDT, a penalty function is used to enforce artificial constraints. Super-convergent patch recovery (SPR) technique and stress-strain relation are used to calculate the in-plane stresses and transverse shear stresses. The FG-GRC material properties are first simplified to isotropic plates and laminated plates, and the results are compared with existing literature to validate the convergence and accuracy of the proposed method. The influences of C⁰-TSDT, C⁰-SSDT, C⁰-ESDT, and the penalty factor on the numerical results are investigated. Furthermore, the effects of graphene nanoplatelets (GPLs) distribution patterns, weight fraction (g_GPL), geometric parameters, etc., on the large deflection bending of FG-GRC plates are discussed.