Nonlocal Strain Gradient Model for the Nonlinear Static Analysis of a Circular/Annular Nanoplate

A nonlinear static analysis of a circular/annular nanoplate on the Winkler–Pasternak elastic foundation based on the nonlocal strain gradient theory is presented in the paper. The governing equations of the graphene plate are derived using first-order shear deformation theory (FSDT) and higher-order shear deformation theory (HSDT) with nonlinear von Karman strains. The article analyses a bilayer circular/annular nanoplate on the Winkler–Pasternak elastic foundation. HSDT while providing a suitable distribution of shear stress along the thickness of the FSDT plate, eliminating the defects of the FSDT and providing good accuracy without using a shear correction factor. To solve the governing equations of the present study, the differential quadratic method (DQM) has been used. Moreover, to validate numerical solutions, the results were compared with the results from other papers. Finally, the effect of the nonlocal coefficient, strain gradient parameter, geometric dimensions, boundary conditions, and foundation elasticity on maximum non-dimensional deflection are investigated. In addition, the deflection results obtained by HSDT have been compared with the results of FSDT, and the importance of using higher-order models has been investigated. From the results, it can be observed that both strain gradient and nonlocal parameters have significant effects on reducing or increasing the dimensionless maximum deflection of the nanoplate. In addition, it is observed that by increasing load values, the importance of considering both strain gradient and nonlocal coefficients in the bending analysis of nanoplates is highlighted. Furthermore, replacing a bilayer nanoplate (considering van der Waals forces between layers) with a single-layer nanoplate (which has the same equivalent thickness as the bilayer nanoplate) is not possible when attempting to obtain exact deflection results, especially when reducing the stiffness of elastic foundations (or in higher bending loads). In addition, the single-layer nanoplate underestimates the deflection results compared to the bilayer nanoplate. Because performing the experiment at the nanoscale is difficult and molecular dynamics simulation is also time-consuming, the potential application of the present study can be expected for the analysis, design, and development of nanoscale devices, such as circular gate transistors, etc.


Introduction
Due to their high stiffness and strength; elastic modulus (higher than 1 TPa); weight ratio; and also outstanding mechanical, electrical, and chemical properties, nanostructures have been applied in diverse domains, including nanoelectromechanical systems, as their high sensitivity and low mass make them perfect for demands in different areas. such as biosensors, medicine, computers, etc. [1][2][3].
The local elasticity model cannot prognosticate extremely small-sized effects on nanostructures because of the deficiency in the nonlocal elasticity model. Because the material's length scale lacks a nonlocal parameter, the classical elasticity model cannot predict the minuscule size role of nanostructures. Experiments and molecular dynamics simulations The vibration behavior of size-dependent structures using the non-local-strain gradient theory or Eringen's theory has been presented in some studies. For example, Qaderi et al. [39] studied cracked graphene-platelet-reinforced composite plates subjected to parametric excitation based on Eringen's theory and the Line spring model. Mahinzare et al. [40] studied the vibration of magnetically actuated viscoelastic functionally graded nanoshells based on the nonlocal strain gradient theory and FSDT. Rashidpour et al. [41] perused the dynamic analysis of the viscoelastic laminated composite nanoplate using the nonlocal strain gradient theory and FSDT, applying the Galerkin method. Ghorbani et al. [42] studied surface effects on the natural frequency of a functionally graded cylindrical nanoshell using the nonlocal strain gradient theory via generalized DQM.
DQM is a powerful numerical method that has been used in many papers [43,44]. The basic goal of this technique is to apply Lagrange interpolation polynomials to field coefficients and to solve the equations at discrete grid points. Improved accuracies can be achieved by using more grid points. Han et al. [45] employed DQM to examine oneelectrode micro-resonators using a generalized 1-DOF model. Liu et al. [46] used DQM to analyze the bending of FGM sandwich plates by considering the tunable auxetic core.
To our best knowledge, FSDT and HSDT with the nonlocal strain gradient theory for the static analysis of circular/annular nanoplates via DQM have not been employed. Furthermore, this paper analyzed a bilayer circular/annular nanoplate located on elastic foundations. The effects of the non-local coefficient, strain gradient, geometric dimensions, boundary conditions, and foundation elasticity on the results of the maximum nondimensional deflection are investigated. Moreover, the deflection results obtained by HSDT have been compared with the results of FSDT, and the importance of using higher-order models has been highlighted. From the results, it can be concluded that both the strain gradient parameter and the nonlocal coefficient have notable effects with respect to reducing (or enhancing) the dimensionless maximum deflection. Additionally, replacing a bilayer nanoplate (considering van der Waals forces between layers) with a single-layer nanoplate (with the same equivalent thickness) is not possible when obtaining exact deflection results. In addition, the deflection of the single-layer nanoplate is lower than the deflection of a bilayer nanoplate. At higher thickness-to-radius ratios, it is better to use HSDT to obtain more accurate results than FSDT. As performing the experiment is difficult because a nanoscale and molecular dynamics simulation is time-consuming, potential applications of the present study can be expected for the analysis and design of nanoscale devices. The results of this paper can be useful for the development of nanostructured devices, such as circular gate transistors, etc.

The Governing Equations for the Axisymmetric Single-Layer Circular/Annular Nanoplate
In this section, equilibrium equations have been derived by employing the energy method for the nonlinear analysis of single-layer circular/annular nanoplates under bending. To derive the equilibrium equations by the energy method, a displacement field using the HSDT and the nonlinear strain components with von Karman's assumptions have been considered. Figure 1 demonstrates a graphene annular plate on the Winkler (k w ) and Pasternak (k p ) elastic foundations with r i (internal radius) and r o (outer radius). Additionally, Figure 2 shows the schematic of the annular plate under bending loads.
Taking into account the HSDT and adding a special function called g(z), the displacement field will be in the form of the following equations which are in the r, θ, and z directions defined by U, V, and W, respectively.   Taking into account the HSDT and adding a special function called g(z), the displacement field will be in the form of the following equations which are in the r, θ, and z directions defined by U, V, and W, respectively.
u0, v0, and w0 are the displacement components of the midplane in the r, θ, and z directions, respectively. Also, ϕ and ψ are the rotation components around the θ and r axes, respectively. The function g(z) is defined as follows.
f (z) and y* are considered different functions used in various references and are listed in   Taking into account the HSDT and adding a special function called g(z), the displacement field will be in the form of the following equations which are in the r, θ, and z directions defined by U, V, and W, respectively.
u0, v0, and w0 are the displacement components of the midplane in the r, θ, and z directions, respectively. Also, ϕ and ψ are the rotation components around the θ and r axes, respectively. The function g(z) is defined as follows.
f (z) and y* are considered different functions used in various references and are listed in  u 0 , v 0 , and w 0 are the displacement components of the midplane in the r, θ, and z directions, respectively. Also, φ and ψ are the rotation components around the θ and r axes, respectively. The function g(z) is defined as follows.
f (z) and y* are considered different functions used in various references and are listed in
The force and momentum resultants in the nonlocal form (NL) are defined as follows:

Derivation of Equilibrium Equations Based on the Energy Method
The potential energy of the system can be defined as the sum of the strain energy caused by the work of internal forces and the potential energy caused by external forces: Π is the potential energy of the entire system, U is the strain energy of the system, and Ω is the potential energy of external forces. According to the principle of minimum potential energy, for a system in equilibrium, the variation in the potential energy is zero: To write the variations of the strain energy of the system, the integral over the volume of strain energy density should be obtained. The strain energy density is as follows: Therefore, for strain energy variations, we have the following: Therefore, δU = V (σ r δε r +σ θ δε θ + σ rz δγ rz )dV (16) Moreover, where k w and k p illustrate the Winkler and Pasternak elastic foundation coefficients, respectively. Therefore, By setting δΠ to zero, the coefficients of δu 0 , δw 0 , and δφ should be zero, and the Euler-Lagrange equations are obtained as follows. All results are nonlocal. Therefore, they are denoted with superscript NL.
Additionally, the third equilibrium equation can be written as follows: The theory of the nonlocal strain gradient (which can be considered as a combination of the nonlocal stress field and strain gradient model) has been derived by Lim et al. [54] and can be expressed as follows: It should be noted that in Equation (23), µ, C ijkl , and l signify nonlocal, elastic, and strain gradient (or internal material length scale) coefficients, respectively. Furthermore, the stress-strain constitutive equation in the nanoscale can be illustrated as follows [55]: It is noted that in Equation (24), E 1 and E 2 denote the Young's modulus along the 1 and 2 directions. Additionally, v 12 and v 21 are Poisson's ratio in the mentioned directions, and G 13 is the shear modulus.
The nonlocal form is described as follows: The force and moment resultants in the local form are as follows: Additionally, the resultants in terms of displacements are obtained in the following forms: The equilibrium equations for a single-layer axisymmetric annular/circular nanoplate on a Winkler-Pasternak elastic foundation are locally expressed in the form of the following equations:

Equilibrium Equations of the Bilayer Axisymmetric Circular/Annular Nanoplate
Graphene sheets have low bending strengths, and to solve this issue, several layers of graphene can be used. In this way, graphene sheets are placed on top of each other and are connected via weak van der Waals bonds, creating layers of graphene [47]. Figure 3 shows the schematic of the bilayer nanoplate under bending load (q) on the elastic foundations. A circular bilayer plate with radius r o and constant thickness h is considered. k 0 is the van der Waals stiffness coefficient between two layers, and kw and k p are the coefficients of the Winkler and Pasternak elastic foundations, respectively.
( ) 2 2 * * 2 13 13 2 The equilibrium equations for a single-layer axisymmetric annular/circular nanoplate on a Winkler-Pasternak elastic foundation are locally expressed in the form of the following equations:

Equilibrium Equations of the Bilayer Axisymmetric Circular/Annular Nanoplate
Graphene sheets have low bending strengths, and to solve this issue, several layers of graphene can be used. In this way, graphene sheets are placed on top of each other and are connected via weak van der Waals bonds, creating layers of graphene [47]. Figure 3 shows the schematic of the bilayer nanoplate under bending load (q) on the elastic foundations. A circular bilayer plate with radius ro and constant thickness h is considered. k0 is the van der Waals stiffness coefficient between two layers, and kw and kp are the coefficients of the Winkler and Pasternak elastic foundations, respectively.  Equilibrium equations for the bilayer circular/annular nanoplate on the Winkler-Pasternak elastic foundation are obtained almost in the same way as for the single-layer nanoplate. The displacement fields for the axisymmetric bilayer circular/annular nanoplate are as follows: index i = 1 represents the first layer, and i = 2 represents the second layer.
The strain equations are similar to those obtained for the single-layer circular/annular nanoplate. Only the energy equation, when using the minimum potential energy principle to obtain the equilibrium equations and boundary conditions for the upper and lower layers, is different according to the following equations: where the upper layer and the lower layer are numbered 1 and 2, respectively. The equilibrium equations in terms of local stresses for the first layer (i = 1) and the second layer (i = 2) are described as follows.

Boundary Conditions
The boundary conditions of the circular/annular plate are as follows: Simply supported (S): Clamped (C):

Dimensionless Assumptions
Due to the small values in the nanoscale and for the simplicity of calculations, the following relations are used to make the governing equations of the nanoplate dimensionless:

Numerical Solution Method
The differential quadrature technique is one of the high-accuracy numerical methods, and it is derived from the quadratic integration method. In the quadratic integration method, the integral at one point along the domain depends on all points along that direction. The value of dependency is determined by weight coefficients.
In the above equation, w 1 , w 2 , . . . , w n are weight coefficients and f 1 , f 2 , . . . , f n are function values at discrete points. Regarding quadratic integration, Belman et al. [56] suggested that the derivative at one point of the function domain depends on the function values at all points of the domain by weight coefficients: In the above equation, A ij is the weight coefficient, and N is the total number of nodes in the direction of r. The weighting coefficients for the first-order derivative are obtained as follows: Additionally, higher-order derivatives are obtained as follows: The weighting coefficients for derivatives of the second and higher orders are introduced in the form of the following equations: In this paper, the distribution of grid points based on Chebyshev-Gauss-Lubato points is used, which enhances the convergence speed of the solution and is described in the following form: In the above equation, r i and r o are the points at the beginning and end of the function, respectively.

Results and Discussion
In this section, different factors are investigated to observe how they influence the deflections of the circular/annular nanoplate based on FSDT and HSDT, and the nonlocal strain gradient model is considered via DQM. In addition, to validate the solution method, the present results (in the case of the circular/annular nanoplate subjected to bending loading loads) are compared with the results of references. Figure 4 shows the effect of the number of nodes used in the differential quadratic method on the results of the present study (the maximum dimensionless deflection of the circular/annular nano plate). As observed, after nine nodes, proper convergence is achieved. In this research study, the number of 11 nodes is used to calculate the results. the following form: In the above equation, i r and o r are the points at the beginning and end of the function, respectively.

Results and Discussion
In this section, different factors are investigated to observe how they influence the deflections of the circular/annular nanoplate based on FSDT and HSDT, and the nonlocal strain gradient model is considered via DQM. In addition, to validate the solution method, the present results (in the case of the circular/annular nanoplate subjected to bending loading loads) are compared with the results of references. Figure 4 shows the effect of the number of nodes used in the differential quadratic method on the results of the present study (the maximum dimensionless deflection of the circular/annular nano plate). As observed, after nine nodes, proper convergence is achieved. In this research study, the number of 11 nodes is used to calculate the results. To validate the current research and solution method, the results obtained have been compared with reference [57], and the results can be observed in Figure 5 for clamped (C) and simply supported (S) boundary conditions. In this reference, the FSDT is used, and the following values are considered :  To validate the current research and solution method, the results obtained have been compared with reference [57], and the results can be observed in Figure 5 for clamped (C) and simply supported (S) boundary conditions. In this reference, the FSDT is used, and the following values are considered: Based on Figure 5, it can be observed that the results of the present paper are in good agreement with the reference. Furthermore, this figure illustrates that the flexibility of the nanoplate decreases with an increasing non-local coefficient [57]. Table 2 compares the deflection of the circular nanoplate with refs. [57][58][59][60] by considering the following assumptions: As can be compared, the results of the present study are in good harmony with the results of the references. Tables 3 and 4 compare the dimensionless maximum deflection obtained from FSDT and HSDT by assuming different functions for a circular nanoplate in clamped and simply supported boundary conditions using different thickness-to-radius ratios (α = h/r o ) and non-local coefficients µ and considering the following: k * w = 1.13 × 10 9 ; k * p = 1.13 × 10 9 , E 1 = 1765 GPa, E 2 = 1588 GPa, ν 12 = 0.3, ν 21 = 0.27, h = 0.34 nm, q * = 1 × 10 9 (67) ering the following assumptions: As can be compared, the results of the present study are in good harmony with the results of the references.   Tables 3 and 4 compare the dimensionless maximum deflection obtained from FSDT and HSDT by assuming different functions for a circular nanoplate in clamped and simply supported boundary conditions using different thickness-to-radius ratios (α = h/ro) and non-local coefficients µ and considering the following: It can be observed in Tables 3 and 4 that by increasing the nonlocal elasticity parameter, the nondimensional deflection is reduced. Moreover, it can be discerned that, by increasing the thickness-to-radius ratio, deflection decreases. It can be observed that for different nonlocal coefficients and thickness-to-radius ratios, the FSDT overestimates the results in comparison with the HSDT.
In this comparison, several different functions have been used as a function to distribute the shear stress along the thickness. These functions are as follows:

Present study (SS) Reference (SS) Present study (CS) Reference (CS) Present study (CC)
Reference (CC) Figure 5. Comparison of the present numerical results with a reference for the maximum dimensionless deflection of the annular nanoplate versus the nonlocal coefficient.   It can be observed in Tables 3 and 4 that by increasing the nonlocal elasticity parameter, the nondimensional deflection is reduced. Moreover, it can be discerned that, by increasing the thickness-to-radius ratio, deflection decreases. It can be observed that for different nonlocal coefficients and thickness-to-radius ratios, the FSDT overestimates the results in comparison with the HSDT.
In this comparison, several different functions have been used as a function to distribute the shear stress along the thickness. These functions are as follows: As observed in Tables 3 and 4, the results of using these functions in the HSDT have very little difference when compared to each other, and it can be concluded that using these various functions interchangeably does not result in a significant difference when calculating nanoscale bending. In fact, compared to other parameters (such as the effect of the nonlocal parameter or strain gradient coefficient), using various g(z) interchangeably has insignificant effects on the results. These functions assume the distribution of the shear stress along the thickness of the plate and satisfy the boundary condition of the shear stress at the top and bottom of the plate. Additionally, unlike the FSDT, there is no need to use a shear correction factor. From the tables, it is observed that deflection decreases with an increase in the order of the theory. At small thicknesses, the difference between the maximum dimensionless deflection of the FSDT and HSDT is almost negligible. As the ratio of the thickness to radius (h/r) increases, the difference between the results of the two theories increases. This significant difference indicates that, on thick plates, more accurate results are obtained by using HSDT. On the other hand, as the nonlocal coefficient is enhanced, the difference between the results of the two theories increases (but this factor is not more significant than the effect of h/r). Figures 6 and 7 illustrate the results of the nondimensional maximum deflection (obtained from HSDT) versus the strain gradient parameter for the circular nanoplate at the clamped and simply supported boundary conditions, respectively. It can be observed that with the enhancement of the strain coefficient, the maximum deflection of the nanoplate decreases. This can be persuaded: by increasing the strain coefficient, the stiffness of the nanoplate is enhanced, so as a result, the deflection of the plate is reduced, which is in agreement with the results in paper [14]. Tables 3 and 4, the results of using these functions in the HSDT have very little difference when compared to each other, and it can be concluded that using these various functions interchangeably does not result in a significant difference when calculating nanoscale bending. In fact, compared to other parameters (such as the effect of the nonlocal parameter or strain gradient coefficient), using various g(z) interchangeably has insignificant effects on the results. These functions assume the distribution of the shear stress along the thickness of the plate and satisfy the boundary condition of the shear stress at the top and bottom of the plate. Additionally, unlike the FSDT, there is no need to use a shear correction factor. From the tables, it is observed that deflection decreases with an increase in the order of the theory. At small thicknesses, the difference between the maximum dimensionless deflection of the FSDT and HSDT is almost negligible. As the ratio of the thickness to radius (h/r) increases, the difference between the results of the two theories increases. This significant difference indicates that, on thick plates, more accurate results are obtained by using HSDT. On the other hand, as the nonlocal coefficient is enhanced, the difference between the results of the two theories increases (but this factor is not more significant than the effect of h/r). Figures 6 and 7 illustrate the results of the nondimensional maximum deflection (obtained from HSDT) versus the strain gradient parameter for the circular nanoplate at the clamped and simply supported boundary conditions, respectively. It can be observed that with the enhancement of the strain coefficient, the maximum deflection of the nanoplate decreases. This can be persuaded: by increasing the strain coefficient, the stiffness of the nanoplate is enhanced, so as a result, the deflection of the plate is reduced, which is in agreement with the results in paper [14].

Non-dimensional deflection
Strian gradient parameter (l) Figure 6. The maximum nondimensional deflection of the circular nanoplate versus the strain gradient parameter (l) in the clamped boundary condition. Table 5 compares the results of FSDT with the average results obtained from HSDT in the clamped boundary condition. It can be seen that by increasing the nonlocal elasticity, the maximum deflection of the nanoplate decreases [57]. Additionally, the results considering different shape functions are almost the same.  Table 5 compares the results of FSDT with the average results obtained from HSDT in the clamped boundary condition. It can be seen that by increasing the nonlocal elasticity, the maximum deflection of the nanoplate decreases [57]. Additionally, the results considering different shape functions are almost the same.   In the diagram, α is defined when thickness h is considered constant (h = 0.34 nm) and the radius changes; moreover, β is assumed when the radius is considered constant (r = 7 nm) and the thickness changes as follows:

Non-dimensional deflection
Strian gradient parameter (l) Figure 7. The maximum nondimensional deflection of the circular nanoplate versus the strain parameter (l) in the simply supported boundary condition. Figure 8 shows the variation of R h (the ratio of the maximum dimensionless deflection obtained from FSDT to HSDT) in terms of α and β ratios for a circular nanoplate with the clamped boundary condition. R h can be defined as follows: where w * FSDT and w * HSDT are the ratios of the maximum dimensionless deflection obtained from FSDT and HSDT, respectively. HSDT) is enhanced with an increasing thickness-to-radius ratio. It can be observed that at the same values of h/r, the differences between these theories are obtained with various α and β ratios. In other words, with the same h/r, the effect of changing the radius on the difference between the results of the two theories is greater than the effect of changing the thickness. Figure 8. Rh versus the thickness-to-radius ratio for α and β ratios in the clamped boundary condition. Figure 9 reveals the nondimensional maximum deflection versus dimensionless bending loads for different nonlocal and strain gradient parameters in the clamped boundary condition. It can be concluded that both the strain gradient and nonlocal parameters have notable effects in reducing or increasing the dimensionless maximum deflection of the nanoplate. In addition, it is observed that with increasing load values, the importance of considering both strain gradient and nonlocal coefficients in the bending analyses of nanoplates is highlighted.  In the diagram, α is defined when thickness h is considered constant (h = 0.34 nm) and the radius changes; moreover, β is assumed when the radius is considered constant (r = 7 nm) and the thickness changes as follows: From the graph, it can be discerned that R h increases as the ratio of the thickness-toradius increases; that is, the difference between the results of the two theories (FSDT and HSDT) is enhanced with an increasing thickness-to-radius ratio. It can be observed that at the same values of h/r, the differences between these theories are obtained with various α and β ratios. In other words, with the same h/r, the effect of changing the radius on the difference between the results of the two theories is greater than the effect of changing the thickness. Figure 9 reveals the nondimensional maximum deflection versus dimensionless bending loads for different nonlocal and strain gradient parameters in the clamped boundary condition. It can be concluded that both the strain gradient and nonlocal parameters have notable effects in reducing or increasing the dimensionless maximum deflection of the nanoplate. In addition, it is observed that with increasing load values, the importance of considering both strain gradient and nonlocal coefficients in the bending analyses of nanoplates is highlighted. Figure 8. Rh versus the thickness-to-radius ratio for α and β ratios in the clamped boundary condition. Figure 9 reveals the nondimensional maximum deflection versus dimensionless bending loads for different nonlocal and strain gradient parameters in the clamped boundary condition. It can be concluded that both the strain gradient and nonlocal parameters have notable effects in reducing or increasing the dimensionless maximum deflection of the nanoplate. In addition, it is observed that with increasing load values, the importance of considering both strain gradient and nonlocal coefficients in the bending analyses of nanoplates is highlighted.

Non-dimensional maximum deflection
Non-dimensional bending load µ=l=0.5 µ=0, l=0.5 µ=0.5, l=0 Additionally, Figures 10 and 11 depict the nondimensional deflection of the annular nanoplate by considering certain nonlocal elasticity and strain gradient coefficients at simply supported (S), clamped (C), and free (F) boundary conditions and by considering the following assumptions: It can be observed that an increase in the nondimensional radius results in a reduction in the nondimensional deflection at all boundary conditions for both layers. Additionally, it can be observed that in the simply supported condition, the deflection is more significantly reduced than in the clamped condition. In other words, the slope of the curve (of the deflection) decreases with decreasing plate flexibility.  It can be observed that an increase in the nondimensional radius results in a reduction in the nondimensional deflection at all boundary conditions for both layers. Additionally, it can be observed that in the simply supported condition, the deflection is more significantly reduced than in the clamped condition. In other words, the slope of the curve (of the deflection) decreases with decreasing plate flexibility.  It can be observed that an increase in the nondimensional radius results in a reduction in the nondimensional deflection at all boundary conditions for both layers. Additionally, it can be observed that in the simply supported condition, the deflection is more significantly reduced than in the clamped condition. In other words, the slope of the curve (of the deflection) decreases with decreasing plate flexibility. Figure 12 shows the variation of the maximum dimensionless deflection in terms of the Winkler and Pasternak elastic coefficients for the circular bilayer nanoplate (the thickness of each layer is 0.34 nm) and the single-layer nanoplate (with a thickness of 0.68 nm) for the clamped boundary condition by considering r = 5 nm, and q = 1 GPa.
forces that hold them together, so the results of a monolayer plate are not the same as the bilayer plate, especially under higher bending loads.
Additionally, it can be observed that, by increasing the stiffness of the elastic foundations, the results converge. In other words, increasing the stiffness of the elastic foundation (or reducing the bending loads) can lead to more similarities between the results obtained from the single-layer plate and the bilayer plate.

Conclusions
The nonlinear bending of the axisymmetric circular/annular nanoplate using the nonlocal strain gradient model with FSDT and HSDT is studied. Additionally, a bilayer circular/annular nanoplate is analyzed. To solve governing equations, DQM was used for the circular/annular nanoplate. For validation, the results were compared with other references. Some of the important results of this paper are as follows:

*
Deflection decreases with the increasing order of the theory. In other words, FSDT overestimates the deflection results, especially at higher thickness-to-radius ratios, so it is better to use HSDT models in those cases. * The ratio of the maximum dimensionless deflection obtained from FSDT to HSDT is enhanced as the ratio of the thickness-to-radius increases. * With the increase in the strain gradient parameter, the nondimensional maximum deflection of the nanoplate decreases. In fact, this figure examines the possibility of replacing a bilayer nanoplate with a single-layer nanoplate. It is clear from the diagram that it is not possible to replace a bilayer nanoplate with a single-layer nanoplate (with the same thickness as the bilayer). In addition, the deflection of the single-layer nanoplate is lower than the deflection of a bilayer nanoplate. In this paper, the interaction between two monolayers of graphene plates is considered a result of the van der Waals interaction. The interpretation (of Figure 12) can be justified as the van der Waals interaction is a relatively weak force. If the bending loads (applied to the plate) are high enough, the molecules break free of the van der Waals forces that hold them together, so the results of a monolayer plate are not the same as the bilayer plate, especially under higher bending loads.
Additionally, it can be observed that, by increasing the stiffness of the elastic foundations, the results converge. In other words, increasing the stiffness of the elastic foundation (or reducing the bending loads) can lead to more similarities between the results obtained from the single-layer plate and the bilayer plate.

Conclusions
The nonlinear bending of the axisymmetric circular/annular nanoplate using the nonlocal strain gradient model with FSDT and HSDT is studied. Additionally, a bilayer circular/annular nanoplate is analyzed. To solve governing equations, DQM was used for the circular/annular nanoplate. For validation, the results were compared with other references. Some of the important results of this paper are as follows: * Deflection decreases with the increasing order of the theory. In other words, FSDT overestimates the deflection results, especially at higher thickness-to-radius ratios, so it is better to use HSDT models in those cases. * The ratio of the maximum dimensionless deflection obtained from FSDT to HSDT is enhanced as the ratio of the thickness-to-radius increases. * With the increase in the strain gradient parameter, the nondimensional maximum deflection of the nanoplate decreases. * Replacing a bilayer nanoplate (considering van der Waals forces between layers) with a single-layer nanoplate (which has the same equivalent thickness as the bilayer plate) is not possible when obtaining exact deflection results, especially when reducing the stiffness of the elastic foundations (or in higher bending loads). In addition, the single-layer nanoplate underestimates the deflection results of the bilayer nanoplate (with the same equivalent thickness as the single-layer nanoplate).
Author Contributions: M.S.: Conceptualization, methodology, software, validation, formal analysis, investigation, and writing-original draft preparation. A.P.: Supervision, conceptualization, methodology, and formal analysis. G.J.: validation and investigation. All authors have read and agreed to the published version of the manuscript.