We combine the advantages of both FGM and GPL using a functionally graded graphene platelet-reinforced composite (FG-GPL RC) structure. GPLs and isotropic epoxy are considered the reinforcement and the matrix of the composite media, respectively. An FG-GPL RC plate with length a, width b, total thickness h, and a circular cutout with radius r, as shown in Figure 1, is assumed in this paper. The plate is resting on a Winkler elastic foundation, where k _w indicates the stiffness of the Winkler elastic foundation. The plate is constructed of N L layers with identical thickness t. GPLs in each individual layer are uniformly distributed, but their volume fraction changes from layer to layer, which leads to a functionally graded pattern. Therefore, each individual layer of the GPL RC plate is isotropic and homogeneous. Figure 2 shows the five uniform or functionally graded distribution patterns of GPL nanofillers. UD denotes the uniform distribution, and FG-X, FG-O, FG-V, and FG-Λ show the functionally graded patterns of GPLs. The UD distribution 1) denotes a homogeneous and isotropic multilayer plate in which nanofillers are uniformly dispersed across the thickness of the plate. The GPL volume fraction is dispersed symmetrically and linearly in FG-X 2) and FG-O 3) distributions, in which the GPL amount gradually changes from the top or bottom surfaces to the middle layer. The unsymmetrical patterns are FG-V 4) and FG-Λ 5) distributions. An adequately large number of layers, N L , leads to a gentle and connected distribution pattern along the plate thickness. The volume fraction V G P L of the kth layer for different GPL patterns is given as follows (Jafari and Kiani, 2021): U D : V G P L k = V G P L * , F G − X : V G P L k = 2 V G P L * 2 k − N L − 1 N L , F G − O : V G P L k = 2 V G P L * 1 − 2 k − N L − 1 N L , F G − V : V G P L k = V G P L * 2 k − 1 N L , (1)

where N L is total number of layers and assumed to be an even number. V G P L * is the average magnitude of the volume fraction of GPLs and can be evaluated based on the weight fraction of GPLs W G P L , density of the matrix ρ m , and GPLs ρ G P L . V G P L * = W G P L W G P L + ρ G P L ρ m − 1 − W G P L . (2)

In the aforementioned equations, the subscripts (m) and (GPL) stand for the matrix and GPLs, respectively. The effective elastic modulus of the nanocomposite plate is obtained through the Halpin–Tsai model, which is dependent on the geometry of the nanoplatelets (Ma and Jin, 2022; Mollaei et al., 2023). Hence, the elasticity modulus of the kth layer of the plate, which is denoted by E k , is regarded as E k = 3 8 1 + ξ L G P L η L G P L V G P L k 1 − η L G P L V G P L k E m + 5 8 1 + ξ W G P L η W G P L V G P L k 1 − η W G P L V G P L k E m . (3) In Eq. 3, the elasticity modulus of the matrix is shown by E m . The two parameters η L G P L and η W G P L in Eq. 3 depend on the elastic properties of the constituents and the geometry of the platelets: η L G P L = E G P L − E m E G P L + ξ L G P L E m , (4) η W G P L = E G P L − E m E G P L + ξ W G P L E m . (5)

The other two parameters, ξ L G P L and ξ W G P L , in Eqs 4, 5 depend on the geometrical dimensions of the nanofillers: ξ L G P L = a G P L t G P L , (6) ξ W G P L = b G P L t G P L , (7) where a G P L , b G P L , and t G P L are assumed to be the length, width, and thickness of GPLs, respectively.

The simple rule of mixture estimation is used to evaluate the Poisson’s ratio and the mass density of the kth layer of the nanocomposite plate (Cong and Duc, 2018; Esmaeili et al., 2022; Farsadi et al., 2022): v k = v G P L V G P L k + v m V m k , (8) ρ k = ρ G P L V G P L k + ρ m V m k , (9)

where V m k = 1 − V G P L k , (10) where v G P L , ρ G P L and v m , ρ m are the Poisson’s ratio and mass density of the GPL and matrix, respectively.

In addition, the rigidity modulus G k of the kth layer is as follows: G k = E k 2 1 + v k . (11)

The governing equations of the plate resting on an elastic foundation are derived applying the third-order shear deformation plate theory. Hence, the displacement field u , v , w is assumed as u x , y , z = u 0 x , y + 5 4 z − 4 3 h 2 z 3 ∅ x x , y + 1 4 z − 5 3 h 2 z 3 w 0 , x , (12) v x , y , z = v 0 x , y + 5 4 z − 4 3 h 2 z 3 ∅ y x , y + 1 4 z − 5 3 h 2 z 3 w 0 , y , (13) w x , y , z = w 0 x , y , (14) where u 0 , v 0 , and w 0 represent the displacements at z = 0 (the midplane of a plate); ∅ x and ∅ y are the transverse normal rotations of the y and x axes. The finite element procedure is used to estimate the eigenvalue buckling forces of the FG-GPL RC plate. A four-node element with seven DOFs per node is employed for modeling the plate. The displacement field of an individual node i is as follows: Q e i = u 0 i v 0 i w i ∅ x i ∅ y i ∂ w 0 i ∂ x ∂ w 0 i ∂ y T , i = 1 , … , 4 . (15)

Among the displacement components, w is estimated by the Hermitian shape functions and the other components by the Lagrangian functions. Hence, the displacement components may be presented as follows: u 0 , v 0 , ∅ x , ∅ y = ∑ i = 1 4 ψ i u o i , ∑ i = 1 4 ψ i v o i , ∑ i = 1 4 ψ i ∅ x i , ∑ i = 1 4 ψ i ∅ y i , (16) w = Ξ 1 w 01 + Ξ 2 ∂ w 01 ∂ x + Ξ 3 ∂ w 01 ∂ y + … + Ξ 10 w 04 + Ξ 11 ∂ w 04 ∂ x + Ξ 12 ∂ w 04 ∂ y , (17) ∂ w ∂ x = ∂ ∂ x Ξ 1 w 01 + Ξ 2 ∂ w 01 ∂ x + Ξ 3 ∂ w 01 ∂ y + … + Ξ 10 w 04 + Ξ 11 ∂ w 04 ∂ x + Ξ 12 ∂ w 04 ∂ y , (18) ∂ w ∂ y = ∂ ∂ y Ξ 1 w 01 + Ξ 2 ∂ w 01 ∂ x + Ξ 3 ∂ w 01 ∂ y + … + Ξ 10 w 04 + Ξ 11 ∂ w 04 ∂ x + Ξ 12 ∂ w 04 ∂ y , (19) where ψ i is the Lagrangian shape function and Ξ i is the Hermitian function. The components of ψ i and Ξ i are presented in Appendix A.

The displacement field of an element is approximated by using its nodal values: U 0 = Δ Η Q e , (20) where Δ Η is the shape function matrix; U 0 and Q e are expressed as follows: U 0 = u 0 , v 0 , w 0 , ∅ x , ∅ y , w 0 , x , w 0 , y T , (21) Q e = Q e 1 Q e 2 Q e 3 Q e 4 T . (22)

The strains of the plate when it is subjected to in-plane pre-buckling loads may be expressed as ε γ = ε 0 γ 0 + z ε 1 0 + z 2 0 γ 2 + z 3 ε 3 0 + ε P G , (23) with ε 0 = ∂ u 0 ∂ x ∂ v 0 ∂ y ∂ u 0 ∂ y + ∂ v 0 ∂ x , ε 1 = 1 4 5 ∂ ∅ x ∂ x + ∂ 2 w ∂ x 2 5 ∂ ∅ y ∂ y + ∂ 2 w ∂ y 2 ∂ ∅ x ∂ y + 2 ∂ 2 w ∂ x ∂ y + ∂ ∅ y ∂ x , ε 3 = − 5 3 h 2 ∂ ∅ x ∂ x + ∂ 2 w ∂ x 2 ∂ ∅ y ∂ y + ∂ 2 w ∂ y 2 ∂ ∅ x ∂ y + 2 ∂ 2 w ∂ x ∂ y + ∂ ∅ y ∂ x , (24) γ 0 = 1 4 5 ∅ x + ∂ w ∂ x 5 ∅ y + ∂ w ∂ y , γ 2 = − 5 h 2 ∅ x + ∂ w ∂ x ∅ y + ∂ w ∂ y ; ε P G = 1 2 ∂ w ∂ x 2 + z 2 2 ∂ ∅ x ∂ x 2 + z 2 2 ∂ ∅ y ∂ x 2 1 2 ∂ w ∂ y 2 + z 2 2 ∂ ∅ x ∂ y 2 + z 2 2 ∂ ∅ y ∂ y 2 ∂ w ∂ x ∂ w ∂ y + z 2 ∂ ∅ x ∂ x ∂ ∅ x ∂ y + z 2 ∂ ∅ y ∂ x ∂ ∅ y ∂ y 0 0 . (25)

By substituting Eq. 15 into Eqs 24, 25, the strain components may be presented as ε = Δ 1 + Δ 2 + Δ 3 Q e , γ = Δ 4 + Δ 5 Q e , (26) with Δ 1 = ∑ i 4 ψ i , x 0 0 0 0 0 ψ i , y 0 0 0 ψ i , y ψ i , x 0 0 0 , (27) Δ 2 = 1 4 ∑ i = 1 4 0 0 Ξ 3 i − 2 , x , x 5 ψ i , x 0 Ξ 3 i − 1 , x , x Ξ 3 i , x , x 0 0 Ξ 3 i − 2 , y , y 0 5 ψ i , y Ξ 3 i − 1 , y , y Ξ 3 i , y , y 0 0 2 Ξ 3 i − 2 , x , y 5 ψ i , y 5 ψ i , x 2 Ξ 3 i − 1 , x , y 2 Ξ 3 i , x , y , (28) Δ 4 = 5 4 ∑ i = 1 4 0 0 Ξ 3 i − 2 , x ψ i 0 Ξ 3 i − 1 , x Ξ 3 i , x 0 0 Ξ 3 i − 2 , y 0 ψ i Ξ 3 i − 1 , y Ξ 3 i , y , (29) Δ 5 = − 5 h 2 ∑ i = 1 4 0 0 Ξ 3 i − 2 , x ψ i 0 Ξ 3 i − 1 , x Ξ 3 i , x 0 0 Ξ 3 i − 2 , y 0 ψ i Ξ 3 i − 1 , y Ξ 3 i , y . (30)

The stress–strain relationships based on Hooke’s law are as follows: σ = D m ε 0 + z ε 1 + z 3 ε 3 , τ = D s γ 0 + z 2 γ 2 , (31) with σ = σ x , σ y , τ x y T , τ = τ x z , τ y z T , (32) D m = E 1 − v 2 1 v 0 v 1 0 0 0 1 − v 2 , D s = E 2 1 + v 1 0 0 1 . (33)

The elastic strain energy of the plate and Winkler elastic foundation is written as follows: U = 1 2 ∑ N E ∫ V e ε T σ   d V + 1 2 k w ∑ N E ∫ S e w 2 d S (34)

or in the matrix form U = 1 2 ∑ N E Q e T K e p Q e + 1 2 ∑ N E Q e T K e w Q e , (35) where the stiffness matrices corresponding to the strain energy of deformation of element K e p and elastic foundation K e w are given as K e p = ∫ S e Δ 1 T A Δ 1 + Δ 1 T Β Δ 2 + Δ 1 T E Δ 3 + Δ 2 T Β Δ 1 + Δ 2 T D Δ 2 + Δ 2 T F Δ 3 + Δ 3 T E Δ 1 + Δ 3 T F Δ 2 + Δ 3 T F Δ 2 + Δ 3 T H Δ 3 + Δ 4 T A ′ Δ 4 + Δ 4 T Β ′ Δ 5 + Δ 5 T Β ′ Δ 4 + Δ 5 T D ′ Δ 5 d S , (36) K e w = ∫ S e Δ 6 T k w Δ 6 d S , (37) Δ 6 = 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 0 Ξ 1 0 0 Ξ 2 Ξ 3 . . . 0 0 Ξ 10 0 0 Ξ 11 Ξ 12 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 . . . 0 0 0 0 0 0 0 0 0 Ξ 1 , x 0 0 Ξ 2 , x Ξ 3 , x . . . 0 0 Ξ 10 , x 0 0 Ξ 11 , x Ξ 12 , x 0 0 Ξ 1 , y 0 0 Ξ 2 , y Ξ 3 , y . . . 0 0 Ξ 10 , y 0 0 Ξ 11 , y Ξ 12 , y 7 * 24 , (38)

where A , B , D , E , F , H = ∫ − h / 2 h / 2 1 , z , z 2 , z 3 , z 4 , z 6 D m d z , (39) A ′ , B ′ , D ′ = ∫ − h / 2 h / 2 1 , z 2 , z 4 D s d z . (40)

The geometric strain energy due to in-plane pre-buckling loads is obtained as U G = 1 2 ∑ N E ∫ V e σ ^ 0 T ε P G d V . (41)

By substituting the geometric strain into Eq. 41, we obtain U G = 1 2 ∑ N E ∫ S e ε ¯ P G T σ 0 ε ¯ P G d S , (42) where ε ¯ P G = w 0 , x w 0 , y ∅ x , x ∅ x , y ∅ y , x ∅ y , y = 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 ∂ ∂ x 0 0 0 0 0 0 ∂ ∂ y 0 0 0 0 0 0 0 ∂ ∂ x 0 0 0 0 0 0 ∂ ∂ y 0 0 u v w ∅ x ∅ y w 0 , x w 0 , y = Δ P G Q e , (43) σ 0 = d i a g h σ ∼ , h 3 12 σ ∼ , h 3 12 σ ∼ , σ ∼ = σ x 0 τ x y 0 ; τ x y 0 σ y 0 T . (44)

Then, Eq. 42 becomes U G = 1 2 ∑ N E Q e T K e G p Q e , (45)

with K e G p = ∫ S e Δ P G T σ 0 Δ P G d S . (46)

For the buckling problem, after assembly of element matrices, we derive the following equation: K P + K w − λ b K G P Q = 0 , (47) where K P , K w , and K G P are the global stiffness matrix, stiffness matrix of the elastic foundation, and global geometric stiffness matrix, respectively. Eq. 47 should be solved to obtain the buckling mode shapes and buckling forces λ b .

In order to validate the accuracy of the results of the present study, the critical buckling load ratios N c r C / N c r M are recalculated and compared as shown in Table 2 with those from Song et al. (2017). A comparative study is carried out for the critical buckling load ratios of a simply supported FG-GPL RC plate without cutouts and with the same material properties listed in Table 1. Here, N c r C and N c r M are critical buckling loads of the plates with and without GPLs, respectively. Geometric dimensions of the plate are a = 0.45   m , a / b = 1 , and a / h = 10 . The number of layers, N = 10 , and 1.0 wt% GPLs are used for comparing results. It can be seen that the present results are in good agreement with those obtained by Song et al. (2017).

In this section, several numerical examples are demonstrated to study the effects of different factors on the buckling behavior of an FG-GPL RC plate with a cutout.

The effect of different thickness ratios ( b / h ) and GPL patterns on the buckling response of a multilayered FG-GPL RC plate with a circular cutout under uniaxial compression is shown in Table 3. The cutout radius is assumed to be r / b = 0.2 and the plate aspect ratio is a / b = 1 , while the plate is fully clamped at all edges. The first four buckling loads are reported in Table 3. This table reveals that with the same volume fraction and same width-to-thickness ratio, the FG-X distribution pattern results in the highest buckling values, while the FG-O distribution pattern results in the lowest values. Hence, UD and FG-V distribution patterns result in intermediate values. In addition, results show that when the plate slenderness ratio increases, the buckling loads decrease, as expected. The influence of GPL weight fraction and distribution pattern of nanofillers on the buckling loads of the FG-GPL RC plate with a circular cutout is provided in Table 4. The buckling loads of the pure epoxy plate are tabulated for comparison purposes. As can be seen in the table, by dispersing only a low amount of GPL nanofillers, the buckling loads can remarkably improve. By dispersion of only 0.5 wt% and 1.0 wt% GPLs for the FG-X plate, the buckling load enhanced 2.97 and 4.84 times, respectively. The effects of two different boundary conditions and GPL patterns on the buckling loads of the plate under uniaxial compression are compared in Table 5. This table shows the first four buckling loads of the plate for a / b = 1 , b / h = 20 , and circular cutout radius r / b = 0.2 . As can be seen in this table, buckling forces are higher for the boundary conditions that restrict more degrees of freedom (DOFs). Values tabulated in Table 6 are provided to indicate the effect of the a / b ratio and GPL pattern on the buckling load results of the FG-GPL RC clamped plate with a circular cutout under uniaxial load. Results denote that by increasing the plate aspect ratio from 1 to 2, the buckling loads enhanced up to 99.6% for the FG-X plate, but by changing the aspect ratio from 2 to 4, the buckling loads decreased by 7.8%. This can be justified by investigating the pre-buckling stresses of the plates with different aspect ratios. In this case, for plates with a/b = 4, the magnitude of the compressive stresses around the hole is more than the plates with a/b = 2. Results of Tables 3–6 denote that the stiffness of the outer layers generally has the most remarkable effect on the stiffness of the composite structure. Therefore, as shown previously, the structure with the FG-X gradient pattern shows the highest stiffness and buckling capabilities, and UD, FG-V, and FG-O patterns in the sequence have lower values.

Table 7 presents the buckling loads for three different radiuses of the circular cutout of the FG-GPL RC plate under uniaxial load. As expected, with the reduction in the cutout radius, the buckling loads are increased. For example, by increasing the cutout radius-to-width ratio ( r / b ) from 0.1 to 0.3, the buckling loads are decreased by 22.3%. This is due to the fact that for larger cutouts, the pre-buckling compressive stresses are more than smaller cutouts. Table 8 lists the first four buckling loads for the FG-GPL RC plate under uniaxial load with a circular cutout resting on the Winkler-type elastic foundation. Values are provided for several elastic foundation stiffness parameters and a plate with CCCC supports. The comparison of the values in Table 8 with those without an elastic foundation (Table 3) denotes that applying an elastic foundation postpones the buckling phenomenon. For example, the buckling loads are increased by 42.6% for the FG-GPL RC plate when the elastic foundation stiffness is increased from 0 to 1 G N m 3 .

Buckling loads of the FG-GPL RC plate under uniaxial and biaxial normal loads are compared with those resting on an elastic foundation in Table 9 for two different plate aspect ratios and four GPL patterns. As can be seen in the table, for the plate aspect ratio of 1, the buckling loads of uniaxial and biaxial compressions are almost identical for the plate without an elastic foundation. For the higher plate aspect ratios ( a / b = 2 ), the buckling load of the plate under biaxial load is considerably lower compared to those under uniaxial load. For all aspect ratios, by adding an elastic foundation to the plate, the buckling loads are increased for both uniaxial and biaxial compressions. Values provided in Table 10 are tabulated with the intent of comparing the critical buckling loads of the FG-GPL RC plate under uniaxial load with different shapes of cutout with/without an elastic foundation. The effect of the elliptical cutout orientation of the plate under uniaxial load is also investigated. The elliptical cutout is considered to have the same area as the circular cutout with r / b = 0.2 . Here, c and d are the semi-major and semi-minor radius of the elliptical cutout, and their radius-to-width ratios are assumed to be 0.3 and 0.133 , respectively. As seen in Table 10, if the major radius of the elliptical cutout is parallel to the direction of the applied load, the critical buckling load of the UD pattern is approximately 2.39 times lower compared to those when the semi-major radius of the elliptical cutout is perpendicular to the applied load. Furthermore, in this case, the buckling loads are lower than the case of the plate with the circular cutout. Results show that the influence of the elastic foundation on the plate with an elliptical cutout is more remarkable compared to the plate with a circular cutout. For example, when the plate with an elliptical cutout is embedded in an elastic foundation, the critical buckling load of the plate with the UD pattern is approximately 1.73 times greater compared with the plate without an elastic foundation, while for the plate with a circular cutout, it is approximately 1.37 times higher. Buckling loads of the FG-GPL RC plate under shear loads are compared with those resting on the elastic foundation, as shown in Table 11, for two different plate aspect ratios, different shapes of cutouts, and four GPL patterns. By increasing the aspect ratio from 1 to 2, the critical shear buckling load of the plate without an elastic foundation and for the UD pattern reduces by 1.83 times, while that of the plate resting on an elastic foundation reduces by 1.55 times. As can be seen from this table, the shear buckling loads of plates with circular and elliptical cutouts are almost identical. Figure 3, Figure 4, and Figure 5 show the first four buckling mode shapes of the full-clamped FG-GPL RC plate discussed in previous numerical tables. Mode shape results indicate that the elastic foundation with higher stiffness generates more slender and local buckling waves.