Buckling and Free Vibration of Plates Resting on Elastic Foundation Using a New Strain Based Finite Element

The dynamic and stability analysis of plate structures resting on elastic foundations has remained prominent in the civil, mechanical, marine and aerospace engineering industries because of its relevance to practical applications. In many cases, the supporting medium is typically modelled using Winkler, Pasternak, or Kerr foundations assumptions, all of which have considerable implications on how both natural frequencies and critical buckling loads can be estimated based upon the foundation’s stiffness parameters. This research effort continues to develop new reliable theoretical and numerical modelling methodologies, which can account for the following factors: transverse shearing effects, material gradation and plate/foundation interactions.

In recent years, interest in theories analyzing plate structures under free vibration and buckling has steadily increased. Several researchers have contributed to the development of Mindlin plate finite elements. Among the early studies was the analysis of free vibration of Mindlin plates by Al Janabi et al. [1]. These authors used a nine-node finite element, with shear strain interpolation within the model to analyze flat square plates subjected to various boundary conditions. Lee [2] proposed a different formulation using a four-node finite element with natural transverse shear strains and improved the calculation of the element’s stiffness and mass matrices. Similarly, Liew et al. [3] proposed a meshfree method of the first-order shear deformation theory (FSDT) using transformation methods to define the plate’s boundary conditions.

Despite these contributions, most Mindlin finite element formulations remain displacement-field based and are often sensitive to shear and distorted meshing, particularly when dealing with thick plates or complex boundary conditions. Exact solutions are also generally limited to specific geometries and boundary conditions.

Alongside these developments, strain-based formulations have emerged as a promising alternative for improving the accuracy and efficiency of the finite element method. Early strain-based formulations dealt with curved elements [4]. Subsequent work focused on the analysis of bending plates [5], [6], [7], [8], [9], [10] and three-dimensional elasticity problems [11], [12], [13]. The strain-based approach has made it possible to reduce parasitic shear stresses in lower-order elements while allowing higher-order displacement fields without adding extra degrees of freedom [14]. Belounar et al. [15] developed a strain-based Mindlin quadrilateral plate element in which curvatures and transverse shear strains are approximated independently. Boussem et al. [16] successfully applied this method to study the vibrational response of rectangular plates coupled to a fluid.

However, the application of strain-based finite elements to plates resting on elastic foundations remains largely unexplored.

Much research has been conducted on plates resting on elastic foundations. Using classical plate theory, Lam et al. [17] studied the elastic bending, buckling, and vibration of isotropic rectangular plates on elastic foundations. Matsunaga [18] developed a higher-order theory for thick rectangular plates on elastic foundations. Bahmyari et al. [19] employed a meshless Galerkin approach for plates on two-parameter foundations, while Nguyen et al. [20] investigated the vibration of plates on Pasternak foundations using FSDT.

Nevertheless, most existing studies on vibration analysis are based on displacement-based or other formulation methods, with limited use of strain-based finite elements.

The buckling behavior of isotropic plates on elastic foundations has also received considerable attention. Hosseini-Hashemi et al. [21], Shufrin and Eisenberger [22], Mizusawa [23], Liew et al. [24], and Akhavan et al. [25] provided valuable analytical and numerical solutions for Mindlin plates under various boundary conditions and foundation models.

However, studies combining buckling analysis of plates on Winkler, Pasternak, or Kerr foundations with strain-based finite element formulations remain very rare.

Despite significant advances in research, this field still faces numerous limitations. Many analytical models are restricted to specific geometries or boundary conditions, while a large proportion of displacement-based finite elements remain highly sensitive to shear and mesh distortion. Furthermore, although higher-order finite elements can improve accuracy, they often require unnecessary additional degrees of freedom or internal nodes, resulting in a substantial increase in computational cost.

To overcome these drawbacks, the present work adopts a strain-based formulation. This approach allows the development of higher-order elements without resorting to internal nodes or additional degrees of freedom, resulting in simple, efficient, and numerically stable elements. In this context, a new four-node quadrilateral strain-based element (SBQ12) is developed within the framework of Reissner-Mindlin plate theory. This element has only three essential degrees of freedom per node (the transverse displacement w and the two rotations $\beta_x$ and $\beta_y$). It relies on an independent interpolation of bending and transverse shear strains, which significantly reduces shear sensitivity for both thin and thick plates.

Furthermore, although several studies have addressed plates resting on elastic foundations, the literature reveals a significant gap regarding the use of strain-based finite element approaches for the free vibration and linear buckling analysis of plates on Winkler, Pasternak, and Kerr foundations. Addressing this gap constitutes the main novelty of this contribution. The free vibration and buckling behavior of plates resting on different elastic foundations are studied in a complete way using a strain-based method. The performance of the proposed SBQ12 element is verified by comparing its results with the existing analytical and numerical solutions using several numerical examples.

The proposed SBQ12 element is a reliable tool that can effectively support experimental research, even though this study is purely numerical. Its predictions of natural frequencies and critical buckling loads on elastic foundations allow the optimization of experimental design, the interpretation of measured structural responses and the calibration of numerical models using experimental data.

In the finite element method, the continuous plate domain is discretized into smaller finite elements. The displacement field within each element is approximated using interpolation (shape) functions that are expressed in terms of the nodal displacements.

The strain-displacement relations for the Reissner/Mindlin plate element (Figure 1) are given by:

where, $\kappa_x, \kappa_y$, and $\kappa_{x y}$ are the curvature strains and $\gamma_{y z}, \gamma_{x z}$ are the transverse shear strains.

In matrix form, it can be written as:

The strain components cannot be considered independent; they must satisfy the compatibility equations, which are given as:

The displacement field resulting from the three rigid body modes is determined by setting Eq. (1) to zero, yielding the following results:

where, $\alpha_1, \alpha_2$ and $\alpha_3$ are three constants.

The proposed quadrilateral element SBQ12 has three degrees of freedom ($w, \beta_x$, and $\beta_y$) at each of the four nodes. Therefore, the displacement field should contain twelve independent constants, and having used three ($\alpha_1, \alpha_2, \alpha_3$) for the representation of the rigid body modes, the remaining nine constants ($\alpha_4, \alpha_5, \ldots, \alpha_{12}$) are to be apportioned among the five assumed strains of the element. The interpolation of the assumed strain field for the present element SBQ12 is expressed as:

The strain vector of Eq. (5) satisfies the compatibility equation (Eq. (3)).

The strain functions($\kappa_x, \kappa_y, \kappa_{x y}, \gamma_{x z}, \gamma_{y z}$) from Eq. (5) are inserted into Eq. (1), and following integration, the following results are obtained:

The displacement functions derived from Eq. (6) are combined with the rigid body mode displacements from Eq. (4) to produce the final displacement shape functions:

The displacement functions of Eq. (7) can be given in matrix form as:

where, $\left[C_u\right]$ is a $3 \times 12$ matrix. $\{\alpha\}$ is the constant parameters vector containing the constants $\left\{\alpha_i\right\} i=1 \ldots 12$.

The following vector $\left\{q_e\right\}$ of total nodal degrees of freedom (dof) for a generic elementary domain $\Omega_e$ was considered:

The vector $\left\{q_e\right\}$ is expressed as a function of the constant parameters $\{\alpha\}$ using Eq. (8) as:

where, $[C]$ is a $12 \times 12$ matrix calculated from the matrix $\left[C_u\right]$ using the coordinates of each node. The constant parameters vector $\{\alpha\}$ can therefore be derived from $E q$. (11):

Substituting Eq. (12) into Eq. (9) yields:

where,

where, $[N]$ is the interpolation functions.

For Reissner-Mindlin plate theory, the curvatures $\{k\}$ and the transverse shear strains $\{\gamma\}$ are given in terms of displacements as:

In which the bending $\left[Q_b\right]$ and the transverse shear $\left[Q_s\right]$ strain matrices are:

The strain vectors can be redefined using the Eq. (12) as:

The element stiffness matrix $\left[K_e\right]$ is given by:

where, $\left[K_e^b\right]$ is the bending part, $\left[K_e^s\right]$ is the shear part and $\left[K_e^f\right]$ is the elastic foundation part which can be expressed in various ways according to the foundation discussed (Winkler [26]/Pasternak [26]/Kerr foundation [27]) as shown in Figure 2:

The element mass matrix $\left[M_e\right]$ is defined as:

The geometrical strains can be expressed as:

where,

$ \left[C_g\right]=\left[\begin{array}{ccc} \frac{\partial}{\partial x} & 0 & 0 \\ \frac{\partial}{\partial y} & 0 & 0 \\ 0 & \frac{\partial}{\partial x} & 0 \\ 0 & \frac{\partial}{\partial y} & 0 \\ 0 & 0 & \frac{\partial}{\partial x} \\ 0 & 0 & \frac{\partial}{\partial y} \end{array}\right]\left[C_u\right] $

The geometrical strains vector is expressed from the dof vector $\left\{q_e\right\}^T$ as:

The element geometrical matrix $\left[K_{e g}\right]$ is formulated as:

The stress–strain relationship is given by:

where, $[D],\left[D_b\right]$, and $\left[D_s\right]$ are, respectively, rigidity, bending rigidity, shear rigidity matrices and $[T]$ is the matrix containing the material mass density.

$ \begin{aligned} & {[D]=\left[\begin{array}{cc} {\left[D_b\right]} & 0 \\ 0 & {\left[D_s\right]} \end{array}\right], \quad\left[D_b\right]=\frac{E h^3}{12\left(1-v^2\right)}\left[\begin{array}{ccc} 1 & v & 0 \\ v & 1 & 0 \\ 0 & 0 & \frac{(1-v)}{2} \end{array}\right],[D_s]=k h G\left[\begin{array}{ll} 1 & 0 \\ 0 & 1 \end{array}\right],} \\ & {[T]=\rho\left[\begin{array}{ccc} h & 0 & 0 \\ 0 & \frac{h^3}{12} & 0 \\ 0 & 0 & \frac{h^3}{12} \end{array}\right],\left[\sigma_0\right]=\left[\begin{array}{cc} \sigma_x^0 & \sigma_{x y}^0 \\ \sigma_{x y}^0 & \sigma_y^0 \end{array}\right],[\tau]=\left[\begin{array}{ccc} h\left[\sigma_0\right] & 0 & 0 \\ 0 & \frac{h^3}{12}\left[\sigma_0\right] & 0 \\ 0 & 0 & \frac{h^3}{12}\left[\sigma_0\right] \end{array}\right]} \end{aligned} $

where, $\sigma_x^0, \sigma_y^0$ and $\sigma_{x y}^0$ are the in-plane stresses.

The matrices $\left[K_e\right],\left[M_e\right]$ and $\left[K_{e g}\right]$ given in Eqs. (22), (28), and (31) are numerically computed. They are then assembled to obtain the structural stiffness, mass and geometrical matrices ($[K],[M]$, and $\left[K_g\right]$).

For free vibration, we use:

For the buckling analysis, we use:

The proposed SBQ12 element was implemented in a MATLAB environment. All element matrices (stiffness $\left[K_e\right]$, mass $\left[M_e\right]$, geometric $\left[K_G\right]$, and foundation stiffness $\left[K_f\right]$ ) are computed at the element level using a full $2 \times 2$ Gauss quadrature scheme, which was found sufficient to achieve exact integration for the present formulation.

The global stiffness, mass, and geometric matrices are assembled using the standard direct stiffness assembly procedure. Boundary conditions are imposed by removing or penalizing the corresponding degrees of freedom in the global matrices. For plates resting on elastic foundations, the foundation contribution $\left[K_f\right]$ is directly added to the global stiffness matrix $[K]$ at the element level.

Free vibration and linear buckling analysis are performed by solving the generalized eigenvalue problem given in Eqs. (33) and (34).

The calculation procedure can be summarized as follows:

1. Definition of the geometry and properties of materials and foundations.

2. Calculation of the element matrices based on the assumed strain fields

3. Assembly of the global matrices using the standard direct stiffness assembly procedure, ensuring compatibility at nodes shared between adjacent elements.

4. Application of boundary conditions using the direct elimination method.

5. Solving the generalized eigenvalue problem for vibration or buckling using MATLAB’s built-in eig function.

6. Post-processing of natural frequencies, vibration modes, and critical loads.

This implementation guarantees the complete reproducibility of all the numerical results presented in this study.

In order to impose different boundary conditions at ($x=0$) and ($x=b$), the following prescribed conditions are considered:

Clamped edges for ($x=0$) and ($x=b$)

$ w=\beta_x=\beta_y=0 $

Simply supported edges for ($x=0$) and ($x=b$)

$ w=\beta_y=0 $

The boundary conditions of the plate are specified by the letter symbols. For instance, two opposite edges simply supported, one edge clamped, one edge free (SCSF) represents a plate with simply supported (S) boundary condition at ($y=0$) and ($y=a$), clamped ($C$) boundary condition at ($x=0$) and free (F) boundary condition at ($x=b$).

Seven combinations of boundary conditions were considered: all edges simply supported (SSSS), all edges clamped (CCCC), two opposite edges simply supported, two opposite edges clamped (SCSC), two opposite edges simply supported, two opposite edges free (SFSF), three edges simply supported, one edge free (SSSF), SCSF, and three edges clamped, one edge free (CCCF), where S, C and F stand for the simply supported, clamped and free boundary conditions, respectively.

The objective of this section is to verify the accuracy and efficiency of the developed SBQ12 quadrilateral element for analyzing linear free vibration and buckling of square and rectangular isotropic plates. Numerous numerical examples were performed on isotropic plates with varying aspect ratios $(a / b)$, thickness ratios $(a / h)$, and boundary conditions to evaluate the performance of the proposed element. The influence of elastic foundations, such as the Winkler, Pasternak, and Kerr models, was also investigated.

To assess the accuracy and convergence behavior of the proposed SBQ12 element, a mesh convergence study was performed on thin square plates ($h / a$ = 0.005) under SSSS boundary condition. The results are compared with the exact solution [28].

Table 1 demonstrates the excellent convergence of the SBQ12 element. With progressive mesh refinement, the results rapidly converge to the exact analytical solution [28]. On a 20 × 20 mesh, the SBQ12 element achieves high accuracy for all modes. This rapid convergence and minimal error confirm the robustness and absence of shear locking of the proposed strain-based formulation.

The validation of the proposed element is carried out by comparing the obtained results with several wellknown references for different test cases. This example considers an isotropic square plate with different boundary conditions. The dimensionless natural frequencies are evaluated for SSSS and CCCC square plates, with different thickness ratios $a / h$ = 10 and 100 , and the first twelve vibration modes were calculated for both plates. The effects of shear deformation are accounted for, and shear correction factors are applied to compare with results from other models, with the shear correction factor set to $\kappa$ = 0.8601.

A dimensionless frequency parameter is defined as:

where, $\omega_{m n}$ is the frequency, a is the plate side length, $\rho$ is the mass density per unit volume, $G$ is the shear modulus given by $G=E / 2(1+v)$.

The results in Table 2 and Table 3 demonstrate that the SBQ12 element effectively avoids shear locking in thin plates ($a / h$ = 100). The relative errors remain generally below 3.5%. This confirms the excellent performance of the strain-based approach in eliminating parasitic shear stresses for thin configurations.

For thick plates ($a / h$ = 10), the element shows high accuracy and stability with maximum errors of 2.37% (SSSS) and 0.96% (CCCC). The low values of error indicate that the SBQ12 formulation is stable in the presence of significant transverse shear deformation.

Overall, the SBQ12 element has a similar accuracy to other elements: strain-based quadrilateral Mindlin plate (SBQMP), four-node plate element developed using the assumed natural strains (ANS4), and Meshfree, while clearly outperforming the R4 formulation, especially for thin plates.

On the other hand, the obtained results highlight the strong influence of boundary conditions and the thickness ratio $a / h$ on the dimensionless natural frequencies of the plate.

Regarding the boundary conditions, CCCC exhibits considerably higher dimensionless frequencies than SSSS for all modes. This result is expected, as the fixed attachment of all edges stiffens the structure and restricts rotation at the boundaries, thereby stiffening the plate and increasing the frequencies.

Regarding the thickness ratio, the transition from a thick plate ($a / h$ = 10) to a thin plate ($a / h$= 100) results in a significant decrease in $\omega^{-}$ for both boundary conditions. This is physically consistent: as the plate becomes thinner relative to its dimension in the a-plane, its flexural stiffness decreases, which lowers $\omega_{m n}$ and consequently $\varpi$.

Furthermore, Figure 3 shows the first four vibration modes of a square SFSF plate with thickness ratio $a / h$ = 10, obtained using the proposed SBQ12 element with a 24 × 24 mesh. These modes clearly reflect the influence of the two opposing free edges. Mode 1 is a basic symmetric bending mode with one half-wave. Mode 2 shows an antisymmetric shape with large displacement along the free edges. Modes 3 and 4 are more complex, with several half-waves and significant twisting and deformation at the free edges. These are typical behaviors of free-edge plates where absence of rotational and transverse stresses provides more flexibility. The SBQ12 element reproduces these characteristic modes accurately.

• Winkler foundation

After validating the SBQ12 element without foundations, its performance was evaluated in more detail by integrating the elastic foundation models of Winkler, Pasternak and Kerr on square and rectangular plates under various boundary conditions (SSSS, SCSC, CCCC, SFSF) and thickness ratios.

Table 4 presents the dimensionless frequencies of square plates resting on a Winkler foundation for different boundary conditions and thickness ratios (0.05, 0.1, 0.15, and 0.2).

The SBQ12 results are in close agreement with the reference solutions from Li et al. [29], the quasi-3D theory [30], and the third-order shear deformation theory (TSDT) [31], with relative errors remaining below 0.6% in nearly all cases. These findings confirm the proper incorporation of the Winkler foundation stiffness term into the element formulation.

The results show that increasing the Winkler stiffness parameter ($k_w$) from 0 to 100 leads to a consistent rise in the dimensionless frequencies, with an increase of approximately 2.4% to 3.0%. This stiffening effect is slightly more noticeable for thinner plates. In addition, the dimensionless frequencies increase significantly with the thickness ratio ($h/a$). Thicker plates exhibit higher dimensionless frequencies due to their greater bending stiffness. These observations demonstrate the combined influence of the foundation stiffness and the plate thickness on the dynamic behavior of the plate-foundation system.

As expected, boundary conditions play a major role: SCSC plates exhibit frequencies 35 to 45% higher than SSSS plates, while SFSF plates exhibit frequencies 45 to 60% lower than SSSS plates. These results highlight the strong sensitivity of the dynamic response to boundary conditions: the more the edges are constrained, the more the overall stiffness of the plate increases, leading to higher dimensionless frequencies.

• Pasternak foundation

The two-parameter Pasternak model, which considers the normal and shear effects of the foundation, was studied on thin square plates ($h/a$ = 100) (Table 5 and Table 6). Compared to the reference results of the FSDT [20] and the element free Galerkin method (EFG) [19], the SBQ12 element achieves good accuracy, with relative errors in the range of 0.2% to 1.7%. The largest errors are observed for CCCC plates with high foundation stiffness, which are mainly due to the increase of shear stresses at the fixed edges. The results clearly indicate the additional stiffening effect of the shear parameter $k_s$. For instance, for $k_w$=100, increasing $k_s$ from 10 to 50 increases the fundamental frequency by about 46.5% for SSSS plates and 25.3% for CCCC plates.

Physically, the Winkler parameter $k_w$ induces a vertical elastic reaction that increases the overall stiffness of the plate-foundation system and hence the natural frequencies. The Pasternak shear parameter $k_s$ provides an additional shear interaction between adjacent springs, which improves the load transfer and decreases local deformations. The combined effect of these two parameters results in a stiffer and more continuous load bearing medium which leads to a substantial increase in dimensionless frequencies.

It is also observed that plates with CCCC exhibit significantly higher natural frequencies than SSSS plate due to the increased rotational restraint, which improves flexural stiffness and limits transverse displacement.

In addition, Table 7 evaluates the effect of the aspect ratio $h/a$ on the dimensionless frequencies of SSSS plates with a fixed thickness ratio $h/a$ = 0.15.

Increasing this ratio $h/a$ from 0.5 to 2 leads to a substantial rise in dimensionless frequencies, reaching up to +128.4% in the absence of elastic foundation. This increase is attributed to the higher bending stiffness in the longitudinal direction. The presence of the Pasternak foundation reduces the sensitivity of the frequencies to the aspect ratio, with a more moderate increase of 92.6%.

Furthermore, Figure 4 shows the effect of this parameter ($a/b$) on the dimensionless natural frequency for SSSS plates ($h/a$ = 0.15) with and without Pasternak foundation ($k_w$ = 100, $k_s$ = 10). The elastic foundation increases the dimensionless frequencies by 60.4%, 36.2%, and 14.9% for aspect ratios a/b of 0.5, 1, and 2, respectively.

Physically, increasing the aspect ratio $a/b$ increases the flexural stiffness in the longitudinal direction, which explains the overall increase in dimensionless frequencies. However, as the plate elongates, its inherent geometric stiffness becomes predominant, thus reducing the relative contribution of the elastic foundation. This explains why the stiffening effect of the foundation is significantly more pronounced for compact plates (low aspect ratio $a/b$) than for elongated plates.

• Kerr foundation

Table 8 studies the performance of the SBQ12 element according to Kerr’s three-parameter elastic foundation model. The analysis focuses on SSSS with a thickness ratio varying from $h / a$ = 0.05 to 0.2 , considering different combinations of normal stiffness $k_u$ and shear stiffness $k_s$, with the bottom layer parameter fixed at $k_l$ = 100.

The Kerr foundation model is a three-parameter elastic foundation that combines the advantages of the Winkler and Pasternak models. The upper spring stiffness $k_u$ represents the normal (vertical) reaction of the foundation, similar to the Winkler modulus. The shear layer stiffness $k_s$ accounts for the interaction between adjacent springs, as in the Pasternak model. The lower spring stiffness $k_l$ connects the shear layer to a rigid base, thus ensuring greater continuity and depth to the foundation’s response.

This three-parameter formulation offers greater flexibility in modeling real soil behavior, especially for foundations with a stiff stratum at a certain depth. It allows a more realistic representation of the continuous interaction between the plate and the supporting medium compared to the one-parameter (Winkler) or two-parameter (Pasternak) models.

The SBQ12 results are in close agreement with the reference solutions of Cao and Vu [37] and quasi-3D theory [30] for all thickness ratios and foundation parameters. The relative errors are very small, remaining below 0.6% for moderate foundation stiffness ($k_u$ = 100, $k_s$ = 0) and generally below 1.5%, even for higher stiffness values. Accuracy is particularly high for thin plates and tends to decrease slightly for thick plates, which is expected given the increased importance of shear deformation.

These results confirm that the SBQ12 element successfully predicts the behavior of the plate resting on Kerr foundation, demonstrating high reliability and numerical stability even with this more complex foundation representation.

The results obtained indicate that the dimensionless frequencies increase significantly with the plate thickness ratio ($h/a$). This trend reflects the greater flexural stiffness of thicker plates. Regarding the Kerr foundation parameters, both the top spring stiffness ($k_u$) and the shear layer stiffness ($k_s$) contribute to increasing the overall stiffness of the system. Higher values of $k_u$ and $k_s$ lead to a steady increase in dimensionless frequencies, with the most significant increase observed when both parameters are increased simultaneously.

These results highlight the important role of the three-parameter Kerr model for modeling the effect of elastic foundation stiffening on plate vibrations.

To evaluate the accuracy and convergence behavior of the proposed SBQ12 element under buckling, a convergence study was carried out on Mindlin square plates under various boundary conditions (SSSS, SCSF, SSSF, SFSF) and different loading cases (uniaxial in the $x$ direction, uniaxial in the $y$ direction and biaxial compression).

As shown in Table 9, the SBQ12 element exhibits rapid and stable convergence with mesh refinement. The results obtained effectively approach the reference solutions of Mizusawa [23], Hosseini et al. [21], and Liew et al. [24]. With a 16 × 16 mesh, the SBQ12 predictions are already very close to the reference values in most cases, and refinement to a 20 × 20 mesh provides only marginal improvements. The element demonstrates excellent performance for all boundary conditions, including those with free edges, which are generally more difficult to handle.

Buckling instability occurs when a square or rectangular isotropic plate subjected to in-plane compression reaches a critical load level, resulting in a sudden deformation out of its plane. Accurate prediction of this critical load is essential for the safe design of plate structures in various fields.

Table 10, Table 11, and Table 12 present the dimensionless critical buckling loads obtained with the SBQ12 element for square plates under uniaxial ($x$ and $y$ directions) and biaxial compression, considering different thickness ratios ($a / h$ = 5 to 1000) and various boundary conditions.

The dimensionless critical buckling load increases with the thickness ratio $a / h$. This trend is typical of dimensionless design, as a higher $a / h$ ratio corresponds to a thinner plate, with results normalized by the flexural stiffness D. Physically, for a given plate length $a$, increasing the thickness $h$ substantially increases the absolute flexural stiffness and, consequently, the actual critical buckling load.

Furthermore, the type of in-plane loading plays a decisive role. Uniaxial compression (Table 10 and Table 11) leads to considerably higher critical buckling load factors compared to biaxial compression (Table 12). For SSSS plates ($a / h$ = 1000), the critical buckling load factor is approximately 50% lower in biaxial compression than in uniaxial compression. Under biaxial loading, the simultaneous application of compressive stresses in both directions induces earlier instability through a more complex two-dimensional deformation mode. In contrast, uniaxial compression induces buckling primarily in one direction, allowing the plate to sustain higher loads before losing stability.

Additionally, Figure 5 displays the first four buckling mode shapes of a CCCF square plate with thickness ratio $a/h$ = 10, subjected to uniaxial compressive load in the $x$-direction. The results were obtained using the proposed SBQ12 element with a regular 24 × 24 mesh. The out-of-plane displacement wis normalized with respect to its maximum absolute value for each mode to better visualize the deformation patterns.

Mode 1 is the basic buckling mode having a single half-wave in the loading direction. Mode 2 shows an antisymmetric deformation and a significant torsion near the free edge. The higher modes (3 and 4) increase in complexity with several half-waves concentrated at the free edge. These modes clearly show the important role of the free edge on the buckling behavior, with the plate deforming more in the unstressed zone.

Table 13 analyses the effect of the aspect ratio $a/b$ and the thickness ratio $a/h$ on the critical buckling load factor under biaxial in-plane compression for simply supported rectangular plates.

The results of the analysis show that an increase in the thickness ratio ($a/h$) from 10 to 100 leads to a slight increase in the critical dimensionless buckling load, approximately 7.3% for $b/a$ = 1 and 4.7% for $b/a$ = 2. This is a typical trend in dimensionless formulations: thinner plates tend to exhibit higher normalized buckling factors due to the normalization by the flexural stiffness.

However, the aspect ratio ($b/a$) has a much more pronounced effect. The critical dimensionless buckling load decreases significantly as the plate changes from a square ($b/a$ = 1) to a rectangular shape ($b/a$ = 2). Physically, as the plate becomes longer, the buckling wavelength increases in the longitudinal direction, which decreases the overall biaxial compressive strength and results in lower critical loads compared to square plates.

To study the buckling of plates resting on elastic foundations, analyses were carried out on simply supported square and rectangular plates subjected to uniaxial and biaxial compressions, using the foundation models of Pasternak and Kerr (Table 14, Table 15, Table 16, and Table 17).

• Pasternak foundation

Table 14 examines the effect of Pasternak foundation parameters on the critical buckling load factor for rectangular plates ($a/b$ = 0.5 and 1) with thickness ratios $a/h$ = 5, 10 and 100 under uniaxial compression.

The results reveal that the buckling load factor increases with the thickness ratio ($a/h$), indicating that thinner plates exhibit higher normalized buckling resistance. The aspect ratio ($a/b$) also has a significant influence: rectangular plates ($a/b$ = 0.5) achieve considerably higher buckling loads than square plates ($a/b$ = 1) under uniaxial compression. Furthermore, incorporating Pasternak foundation parameters ($k_w$ = 100 and

$k_s$ = 10) significantly improves the buckling capacity compared to the case without foundations. This improvement demonstrates the contribution to stiffening due to both the vertical spring stiffness ($k_w$) and the shear layer stiffness ($k_s$), which together provide additional resistance to buckling deformation.

In the same context, Table 15 presents the influence of the foundation stiffness on the critical buckling load parameter for square plates subjected to uniaxial and biaxial compression under SSSS using the proposed element SBQ12, for a thickness ratio $a/h$ =20. The material properties are: Young’s modulus $E$ = 151 GPa and density $\rho$ = 5700 kg/m$^3$.

The results show that the elastic foundation parameters have a very significant effect on the dimensionless buckling load. The buckling load factor increases by approximately 35% for ($k_w$, $k_s$) = (100, 10), and by 227% for ($k_w$, $k_s$) = (1000, 100). This effect is even more pronounced under biaxial compression, where the buckling load factor can increase by up to 353% for the stiffest foundation configuration ($k_w$, $k_s$) = (1000, 100).

• Parametric study of buckling with Kerr foundation

This section presents a parametric study of the buckling behavior of square and rectangular plates on a Kerr foundation.

This analysis is particularly significant due to the scarcity of available references dealing with buckling analysis of plates resting on Kerr-type foundation.

Table 16 and Table 17 present the dimensionless critical buckling loads for square and rectangular plates ($a/b$ = 0.5 and 1) with different boundary conditions (SSSS, SCSC, and CCCC) under uniaxial and biaxial compression. The analyses cover thickness ratios $a/h$ = 5 to 100, with a fixed lower layer parameter $k_l$ = 100 and various combinations of normal stiffness $k_u$ and shear stiffness $k_s$.

The dimensionless critical buckling load is expressed as:

The results show clearly that increasing the Kerr foundation stiffness parameters ($k_u, k_s$) significantly improves the buckling resistance of both rectangular and square plates. For the rectangular plate ($a / b$= 0.5) and SSSS boundary conditions at a thickness ratio $a / h$ = 5, the dimensionless critical load increases by approximately 83%, confirming the importance of the normal stiffness $k_u$. A similar trend is observed under biaxial loading, where the load increases by approximately 82.8%. For the square plates ($a / b$ = 1), the same comparison reveals an increase of approximately 42% under uniaxial loading and approximately 46% under biaxial loading, indicating a slightly lower, but still significant, sensitivity to foundation stiffness.

Furthermore, the influence of the thickness ratio is evident: increasing $a / h$ from 5 to 100 systematically increases the dimensionless buckling load due to increased flexural stiffness. The influence of boundary conditions is even more pronounced: for ($a / b$ = 0.5, $k_u$ = 100, $k_s$ = 10, $a / h$ = 5), the critical load increases from 2.4413 (SSSS) to 3.9661 (CCCC), representing an increase of approximately 62.5%. Similar trends are observed for other configurations, sometimes exceeding 70%.

From a physical perspective, the increase in $k_u$ strengthens the vertical support, while $k_s$ introduces a shear interaction within the foundation, improving load transfer and reducing local deformations. Combined with greater thickness and higher edge stresses, these effects significantly improve the overall stiffness of the platefoundation system, resulting in better resistance to buckling under uniaxial and biaxial compressive loads.

To further investigate the role of the bottom layer stiffness parameter $k_l$ in the Kerr foundation model, additional calculations were performed for $k_l$ = 50 and $k_l$ = 200 (Table 18, Table 19, Table 20, and Table 21), for square and rectangular plates ($a / b$ = 0.5 and 1) under uniaxial and biaxial compression, maintaining the same plate boundary conditions (SSSS, SCSC, CCCC) and the same combinations of $k_u$ and $k_s$.

The results reveal a clear and consistent influence of $k_l$ on the critical dimensionless buckling load. As $k_l$ increases from 50 to 200, the plate’s buckling resistance generally varies, but the magnitude of this effect is highly dependent on the specific combination of $k_u$ and $k_s$, as well as on the plate geometry (aspect ratio $a / b$) and the type of loading (uniaxial or biaxial).

It is noteworthy that when $k_s$ =10 and $k_u$ = 0, the values of $\widetilde{N}_{c r}$ remain virtually unchanged, regardless of the $k_l$ value. For example, for the SSSS boundary conditions at $a / h=5$, the result is 3.1415 for all three values of $k_l$ (50,100, and 200) for the square plate ($a / b$ = 1). The same invariance is observed for the rectangular plate ($a / b$ = 0.5): for SSSS at $a / h$ = 5, $\widetilde{N}_{c r}$ = 1.3352 regardless of $k_l$. This indicates that when the normal stiffness $k_u$ is absent, the lower layer parameter $k_l$ does not contribute to buckling resistance for either geometry, and only the shear stiffness $k_s$ governs the foundation’s response.

Conversely, when $k_u$ = 100 and $k_s$ = 0, the influence of $k_l$ becomes significant for both plate geometries. In the SSSS case with $k_u$ = 100 and $k_s$ = 0 for $a / h$ = 5, $\widetilde{N}_{c r}$ increases from 3.4673 ($k_l$ = 50) to 3.6301 ($k_l$ = 100), then to 3.7929 ($k_l$= 200), representing a total increase of approximately 9.41% over this range for the square plate ($a / b$ = 1). For the rectangular plate ($a / b$ = 0.5) under the same conditions, $\widetilde{N}_{c r}$ increases from $1.6651\left(k_l=50\right)$ to $1.8300\left(k_l=100\right)$ and further to $1.9950\left(k_l=200\right)$, representing a larger relative increase of approximately $19.81 \%$. This stronger sensitivity in the rectangular case can be attributed to its lower inherent bending stiffness along the longer dimension, which makes the vertical foundation support $k_u$ more influential in resisting out-of-plane deformation.

A similar trend is observed for the SCSC and CCCC boundary conditions of both geometries, although the relative gain decreases with increasing edge restraint, as the plate stiffness is already high in the clamped case. For instance, for the rectangular plate ($a / b$ = 0.5) with CCCC boundary conditions, $k_u$ = 100, $k_s$ = 0, and $a / h$ = 5, $\widetilde{N}_{c r}$ increases from 3.3017 ($k_l$ = 50) to 3.4098 ($k_l$ = 100) and to 3.5149 ($k_l$ = 200), a total gain of only about 6.4%, compared to 19.81% for the SSSS case. This confirms that clamped edges already provide substantial restraint, reducing the marginal contribution of $k_l$.

In the combined case $k_u$ = 100 and $k_s$ = 10, the results show a consistent and monotonic decrease in $\tilde{N}_{c r}$ as $k_l$ increases from 50 to 200 , for all boundary conditions and both aspect ratios. In the SSSS case at $a / h$ = 5, $\widetilde{N}_{c r}$ drops from 4.6227 ($k_l$ = 50) to 4.4666 ($k_l$ = 100) and then to 4.3105 ($k_l$ = 200) for the square plate ($a / b$ =1). The same trend is observed for SCSC and CCCC boundary conditions. For the rectangular plate ($a / b$ = 0.5) under SSSS conditions, $\widetilde{N}_{c r}$ similarly decreases from 2.4802 ($k_l$ = 50) to 2.4413 ($k_l$ = 100) and to 2.4025 ($k_l$ = 200), confirming that this behavior is independent of the plate geometry. Under biaxial loading (Table 19 and Table 21), the same monotonic decreasing trend is observed for both geometries when $k_u$ = 100 and $k_s$ = 10, indicating that the coupled nature of the Kerr parameters is not affected by the type of in-plane loading.

This behavior indicates that when $k_u$ and $k_s$ are active simultaneously, a stiffer lower layer ($k_l$) tends to reduce the overall buckling resistance. This is explained by the mechanical role of $k_l$ in the Kerr model: a higher stiffness of the lower spring reduces the allowable deformation in the lower layer, which limits the engagement and load transmission through the shear layer $k_s$. This effect is present regardless of the plate’s aspect ratio, though its absolute magnitude is smaller for rectangular plates ($a / b$ = 0.5) due to their lower reference buckling loads.

It is also worth noting the influence of the loading type on the sensitivity to $k_l$. Under biaxial compression, the reference buckling loads are roughly half those under uniaxial compression for both plate geometries, as expected from the two-directional stress state. However, the relative trends with respect to $k_l$ remain the same: $k_l$ has no effect when $k_u$ = 0, increases $\widetilde{N}_{c r}$ when only $k_u$ is active, and decreases it when both $k_u$ and $k_s$ are simultaneously active. This consistency across loading types further confirms the robustness of the identified coupling mechanism within the three-parameter Kerr foundation model.

Consequently, the effective contribution of the foundation to the plate stiffness is reduced when $k_l$ is large and $k_s$ is already active, revealing the coupled, rather than additive, nature of the three-parameter Kerr foundation. This coupling effect is consistent across all tested configurations: square and rectangular plates, uniaxial and biaxial compression, and all boundary conditions considered.

This study presents a new strain-based quadrilateral finite element (SBQ12), grounded in Reissner–Mindlin plate theory, for the analysis of isotropic plates resting on elastic foundations. The element is characterized by three degrees of freedom per node (transverse displacement and two rotations) and relies on independent interpolation of bending and transverse shear strains, which naturally avoids the use of internal nodes or additional degrees of freedom.

The proposed formulation, which uses independent interpolations of bending and transverse shear strains, effectively eliminates shear locking and ensures excellent numerical stability for both thin and thick plates.

The numerical results demonstrate the excellent performance of the proposed SBQ12 element, with relative errors generally below 3.5% for free vibration analyses and 4% for buckling analyses across a wide range of boundary conditions, thickness ratios, and foundation stiffnesses. Larger deviations (up to 4%) were occasionally observed for very thick plates or highly restrained boundary conditions, where transverse shear effects become more pronounced.

Overall, the developed SBQ12 element has proven to be accurate, robust, and computationally efficient for the free vibration and linear buckling analysis of isotropic plates resting on Winkler, Pasternak, and Kerr elastic foundations. Its strain-based formulation effectively eliminates shear locking while maintaining a minimal number of degrees of freedom, making it a reliable tool for practical engineering applications.

The parametric studies conducted in this work have highlighted several important physical trends. For the Winkler and Pasternak foundations, increasing the foundation stiffness parameters consistently raises both natural frequencies and critical buckling loads, with the effect being more pronounced for compact plates (low $a / b$ ratio) and simply supported boundary conditions. The Pasternak shear parameter $k_s$ introduces an additional stiffening effect that significantly amplifies the response beyond what the Winkler spring alone can provide: for instance, increasing $k_s$ from 10 to 50 (with $k_w$ = 100) raises the fundamental frequency by up to 46.5% for SSSS plates.

For the three-parameter Kerr foundation, the parametric study over $k_l$ = 50, 100, and 200 revealed a nontrivial coupling mechanism: when only $k_u$ is active ($k_s$ = 0), increasing $k_l$ consistently raises the buckling resistance, with a more pronounced effect for rectangular plates ($a / b$ = 0.5, 19.8%) than for square plates (9.41%). However, when both $k_u$ and $k_s$ are simultaneously active, a stiffer lower layer $k_l$ systematically reduces the critical buckling load, revealing the coupled and non-additive nature of the Kerr foundation parameters. This behavior is consistent across all boundary conditions, plate geometries, and loading types (uniaxial and biaxial compression) considered in this study.

Although the present study provides a comprehensive analysis, some limitations should be acknowledged. The formulation is currently restricted to isotropic plates, linear free vibration, and linear buckling under static in-plane loads. It does not account for nonlinear geometric effects, post-buckling behavior, or material nonlinearity. Furthermore, the results presented are purely numerical, and experimental validation on physical specimens remains to be carried out.

Future research will extend the SBQ12 element to laminated composite plates, functionally graded materials (FGM), and nonlinear problems. Experimental validation on scaled plate-foundation models is also planned to confirm the numerical predictions. Furthermore, the element will be implemented in general-purpose finite element software to facilitate its use by practicing engineers. Extensions to dynamic loading, including transient response and forced vibration under moving loads, are also envisaged as part of future investigations.