Mathematical Models with Buckling and Contact Phenomena for Elastic Plates: A Review

A review of mathematical models for elastic plates with buckling and contact phenomena is provided. The state of the art in this domain is presented. Buckling effects are discussed on an example of a system of nonlinear partial differential equations, describing large deflections of the plate. Unilateral contact problems with buckling, including models for plates, resting on elastic foundations, and contact models for delaminated composite plates, are formulated. Dynamic nonlinear equations for elastic plates, which possess buckling and contact effects are also presented. Most commonly used boundary and initial conditions are set up. The advantages and disadvantages of analytical, semi-analytical, and numerical techniques for the buckling and contact problems are discussed. The corresponding references are given.

plate equations work best for small deflections of the plate ( [23]). They become not so accurate in the case of large displacements of the plate. To describe the large displacements, nonlinear models have been created (e.g., those of Antmann [24], Lagnese [25], and Sokolnikov [26]). The nonlinear equations originate from the Föppl-von Kármán nonlinear model. They have a wide spectrum of applications in applied science and engineering. Solving the nonlinear von Kármán system is difficult in virtue of the high order partial differential equations and presence of nonlinearity. Therefore, a classical solution has not yet been found. However, some semi-analytical and numerical approaches have been developed for solving this problem.
Here, we focus our attention on the plate models, which combine buckling and unilateral contact effects together. The buckling and contact problems are important subjects of investigation in applied mechanics and engineering. Buckling occurs when the plate is compressed on its edges. Most worthwhile is the investigation of postbuckling behavior and bifurcation phenomenon for a solution in the nonlinear equations. Contact effects occur when the plate rests on elastic foundation or when the plate is delaminated. The delamination has a place in the composite plates.

Governing Equations
Let us consider a plate which is under load q in vertical z direction, as shown on Figure 1.
There are two basic classical linear plate models. The first one is based on the Kirchhoff-Love principles and the second one issues from the Mindlin-Reissner theory for thick plates. The linear plate equations work best for small deflections of the plate ( [23]). They become not so accurate in the case of large displacements of the plate. In order to describe the large displacements, nonlinear models have been created (Antmann [24], Lagnese [25], Sokolnikov [26] etc). The nonlinear equations originate from the Föppl -von Kármán nonlinear model. They have a wide spectrum of applications in applied science and engineering. To solve the nonlinear von-Kármán system is difficult in virtue of the high order partial differential equations and presence of nonlinearity. Therefore, it has not been found a classical solution for it up to now. However, some semi-analytical and numerical approaches have been developed for solving this problem.
Here we focus our attention on the plate models, which combine buckling and unilateral contact effects together. The buckling and contact problems are important subjects of investigation in applied mechanics and engineering. The buckling occurs when the plate is compressed on its edges. Most worthwhile is investigation of postbuckling behaviour and bifurcation phenomenon for a solution in the nonlinear equations. Contact effects occur when the plate rests on elastic foundation or when the plate is delaminated. The delamination has a place in the composite plates.

Governing Equations
Let us consider a plate which is under load q in vertical z direction as shown on Figure 1. Below we present the Kirchhoff-Love linear equation and the von-Kármán nonlinear system. Both of them are derived by using the principle of virtual work from the equilibrium equations (Ciarlet [2], Reddy [14], Bloom; Coffin [50] etc). Below, we present the Kirchhoff-Love linear equation and the von Kármán nonlinear system. Both are derived by using the principle of virtual work from the equilibrium equations (e.g., Ciarlet [2], Reddy [14], and Bloom and Coffin [50]).
According to the Kirchhoff-Love linear theory for the deflection w = w(x, y) of two variables x, y with taking into account the vertical loading q, the following equation holds where ∆ 2 is the biharmonic operator and G is the domain, occupied by the plate. In the case of a rectangular plate, G = (0, l 1 ) × (0, l 2 ), where l 1 and l 2 are the lengths of sides of the plate. Further, D = Eh 3 /12(1 − ν 2 ) is the flexural rigidity of the plate, h is the thickness of the plate, E is the Young modulus, and ν is the Poisson coefficient.
In the nonlinear case, the symmetric second rank Green-Saint Venant tensor (the Lagrangian strain tensor) is applied (see, e.g., [9]). The following equation for the stress function ψ = ψ(x, y) holds and for the displacement w taking into account vertical loading q, the following equation reads where [w, ψ] is called the Monge Ampère form, Two fourth-order Equations (2) and (3) are coupled through the nonlinear terms. The considered nonlinear model has been extended to the shallow shell model in the work of Vorovich [51]. Nonlinear functional analysis is used to investigate the stability of the models.
In the linear case, one obtains two uncoupled equations: (1) The Airy equation for plane stress: (2) The linear plate bending Kirchhoff-Love's model in Equation (1) (originated from Germain's equation).

Boundary Conditions
Most often, meet boundary conditions for plates are simply supported, partially clamped, totally clamped, Robin's and free boundary conditions along the edges of the plate (see, e.g., Bloom and Coffin [50], Ciarlet and Rabier [3], and Chien et al. [52]).
If all edges of a plate are simply supported, then the displacement and the Airy stress potential satisfy the following boundary conditions, A more appropriate condition in a physical sense for the stress function is (Schaeffer and Golubitsky [53]) where ∂/∂n denotes the normal derivative. The simply supported boundary conditions are convenient from a mathematical point of view, but they are notoriously difficult to achieve in experiments.
The most realistic are totally clamped boundary conditions, in (0, T] × ∂G for the dynamic case.
The partially clamped boundary conditions are combinations of the simply supported and the clamped boundary conditions. If the clamped boundary conditions are relaxed at the ends of the plate, then such boundary conditions are called the Robin boundary conditions. For example, where g i (µ) ∈ R, i = 0, 1 are smooth functions and satisfy If a plate, for example, has a free edge parallel to the axis y at x = l 1 , then the following boundary conditions are considered

Buckling Phenomenon
In this section, we define buckling and postbuckling effects on example of the nonlinear static equations for the plate. A general analysis of plate buckling and postbuckling behavior is given in the handbook of Bloom and Coffin [50]. In particular, a derivation of von Kármán nonlinear system with buckling is presented.
If, except vertical loads q, in-plane distributed loads are given along the boundary ∂G of the plate as a regular function θ(x, y), defined on G, then Equation (3) for the vertical displacement w yields The values of θ depend on the portion of the boundary, subjected to compression or tension. Various problems can be modeled by θ(x, y). The load can be considered compressive (−) or tensile (+). If we have compressive loads on all the sides of the plate, then where λ is the scalar compressive loading factor (magnitude). Compression, for example, along two x-parallel sides, is denoted If a plate is compressed along two x-parallel sides and the other sides are under tension, then where λ is the tensile loading parameter. Analogous expressions hold for the other cases. The compression and tension of the plate are shown in Figure 2.
If q ≡ 0, then Equation (4) with the most commonly used boundary conditions becomes an eigenvalue problem.
where λ is the tensile loading parameter. Analogous expressions hold for the other cases. The compression and tension of the plate are shown in Figure 2.
If q ≡ 0 then the Equation (4) with the most commonly used boundary conditions become an eigenvalue problem. For the problem (4), (2) the following cases are true four sides by compr. load, ∂ 11 , two x-sides by compr. load, ∂ 22 , two y-sides by compr. load, comp. on x-sides and tens. on y-sides, 11 , comp. on y-sides and tens. on x-sides, where λLw = [θ, w], ∂ 11 = ∂ 2 /∂x 2 , ∂ 22 = ∂ 2 /∂y 2 . Then from the equation (4) and definitions (5) we obtain According to the results of [38] if the compressive loading factor λ is greater than Euler critical value λ E (the first eigenvalue of the linearized problem for the nonlinear equation (6) with the fixed tensile loading (λ = λ f ) then the plate buckles and there exist at least three solutions for (6), (2): (w, ψ), (−w, ψ) for w = 0 and ψ = 0, and (0, 0). Otherwise, the plate is not deformed and it remains flat. The critical buckling loading factor usually depends on Poisson's ratio and Young's modulus. The Euler critical value λ E also depends on the tensile loading factor λ f . If the tensile (stretching) load increases then λ E = λ E (λ ) increases as well. For the problems in Equations (2) and (4), the following cases are true four sides by compr. load, ∂ 11 , two x-sides by compr. load, ∂ 22 , two y-sides by compr. load, comp. on x-sides and tens. on y-sides, 11 , comp. on y-sides and tens. on x-sides, where λLw = [θ, w], ∂ 11 = ∂ 2 /∂x 2 , ∂ 22 = ∂ 2 /∂y 2 . Then, from Equation (4) and the definitions in Equation (5), we obtain According to the results of [38], if the compressive loading factor λ is greater than the Euler critical value λ E (the first eigenvalue of the linearized problem for the nonlinear equation in Equation (6) with the fixed tensile loading (λ = λ f ), then the plate buckles up and there exist at least three solutions to Equations (2) and (6): (w, ψ), (−w, ψ) for w = 0 and ψ = 0, and (0, 0). Otherwise, the plate is not deformed and it remains flat. The critical buckling loading factor usually depends on Poisson's ratio and Young's modulus. The Euler critical value λ E also depends on the tensile loading factor λ f . If the tensile (stretching) load increases, then λ E = λ E (λ ) increases as well.
Below, we discuss the existing semi-analytical and numerical techniques for treating buckling and contact plate problems and the state of the art in this direction.
As it is mentioned above, a classical solution for Equations (2) and (6) has not yet been found. However, the buckling and postbuckling behavior of the solution of the nonlinear model in Equations (2) and (6) has intensively been studied numerically by many authors, e.g., Caloz and Rappaz [42], Dossou and Pierre [49], Chien et al. [44], Chien et al. [52], Muradova [46,47], and Matkowsky and Putnick [43]. The existing techniques for treating the nonlinear mechanical model are mainly based on finite element analysis, finite difference approximations, and spectral and pseudo-spectral methods.
In the works [46,47], the Fourier transform is proposed for Equations (2) and (6). The solution is expanded into double Fourier series and the partial sums are considered. The coefficients in the series are good guides for following branches of the solution and the trigonometrical functions reflect the shape of the eigenfunctions. A detailed bifurcation analysis, based on spectral approach, for the simply supported, partially clamped, and totally clamped plates is presented in [47]. The buckling loads are computed and trained by the neural network for prediction in the paper [54]. A topology of bifurcation diagrams for the von Kármán problem for totally clamped plate is analyzed by finite elements, combined with Newton-GMRES algorithm, in [49] and by the spectral method with numerical continuation in [46]. A finite element method (FEM) together with a predictor-corrector method and block GMRES algorithm was implemented by Chien et al. [44] for tracing curves of the solution for the partially clamped plate. The mode jumping phenomenon, which implies non-stability of the primary solution branches through further bifurcations, for the partially clamped boundary conditions was investigated by Holder and Schaeffer [55], Schaeffer and Golubitsky [53], and later by Chien et al. [45].
In the case of employing spectral or pseudo-spectral approaches, the choice of the global basic functions depends on the type of boundary conditions. For example, for the simply supported, partially clamped, and totally clamped plates, the global basic functions are constructed, based on trigonometrical and Legendre's polynomials. In the case of the presence of free plate edges, a combination of trigonometrical functions and the Chebyshev polynomials can be used. In this case, the spatial discretization can be done through a collocation method.
The spectral and pseudo-spectral methods give a high accuracy of computations while finite elements are more flexible due to the local choice of the bases.

Contact Models with Elastic Foundations
If instead of the transversal loading forces q we have elastic foundation forces in Equation (6), then where p is the reaction of a general subgrade model for nonlinear elastic foundations. The function p is the distributed reaction force per unit area.
The plate comes in contact with the upper or lower foundations (obstacles) if w > 0 or w < 0, respectively (see Figure 3). The plate comes in contact with the upper or lower foundations (obstacles) if w > 0 or w < 0, respectively (see Figure 3).
If p ≡ k 1 w, then we deal with a linear Winkler-type elastic foundation. The coefficient k 1 is the spring constant of the soil. A review of Winkler's foundation and numerous applications of it, including adhesion and historical perspective are given in the paper of Dillard et al [56]. The adaptation and influence of the Winkler's model on the displacement are also discussed. Various extensions of the Winkler's model are given.
The Winkler's spring model usually does not simulate all the properties of the soil. Therefore, two-parameter foundation has been introduced in order to describe interaction among the springs. Most used two parameters foundations are Filonenko-Borodich, If p ≡ k 1 w, then we deal with a linear Winkler-type elastic foundation. The coefficient k 1 is the spring constant of the soil. A review of Winkler's foundation and numerous applications of it, including adhesion and historical perspective, is given in the paper of Dillard et al. [56]. The adaptation and influence of the Winkler's model on the displacement are also discussed. Various extensions of the Winkler's model are given. The Winkler's spring model usually does not simulate all the properties of the soil. Therefore, the two-parameter foundation has been introduced to describe interaction among the springs. The most used two-parameter foundations are Filonenko-Borodich, where T is a constant representing the tension action of a thin elastic membrane, and shear Pasternak-type foundation model where k g is the shear modulus, in others words the interaction parameter. Another type of description is Hetényi model, where S is the flexural rigidity of the soil layer. The foundation models in Equations (7)-(9) and other type models are presented in the works of Dillard et al. [56] and Wang et al. [57]. They are all particular cases of the general subgrade model. For many practical foundations, Equation (10) consists of separated nonlinear elastic Winkler-type foundation springs with contribution k 1 w − k 2 w 3 and shear Pasternak-type foundation contribution k g ∆w. The constants of the foundations depend on the displacement in a non-smooth way.
When contact occurs (p = 0) and the reaction forces of the foundation are evaluated by Equation (10), the Euler critical value is the first eigenvalue of the linearized equation, of the problems in Equations (2) and (4). Various types of contact problems in solid mechanics and numerical treatments of them are considered along with others in the books of Wrigers [58] and Selvadurai [59]; and in the works of Dumir and Bhakar [60], Katsikadelis and Yiotis [61], and Papanikolaou and Doudoumis [62]. The mathematical theory for unilateral constraints including the Kirchoff plates was developed by Fichera [41]. There are few devoted to the investigation of the two phenomena buckling and contact together in the nonlinear cases. Nevertheless, one can mention the articles of Yu and Wang [63], Bielski and Telega [64], Shen [65], Borisovich et al. [34], Muradova et al. [40], and Muradova and Stavroulakis [37,38]. The works by Borisovich et al. [34] and Muradova et al. [40] are devoted to the exploration of stability of a rectangular, isotropic plate on the linear Winkler-type foundation. In [40], an application of the contact model to the problem from biomechanics and the numerical treatment of it are presented. Buckling and postbuckling analysis of the simply supported plate with the nonlinear foundation of the type in Equation (10) is studied in the work [65] and also in [37,38]. In the paper of Holanda and Goncalves [66], the buckling and postbuckling behavior of plates laterally constrained by a tensionless foundation and subjected to in-plane compressive forces are investigated. Mechanical models for buckling of unilaterally constrained by nonlinear elastic foundations rectangular and infinite plates are given in the papers of Shahwan and Waas [67,68], respectively. In [67], the influence of different boundary conditions for the rectangular plates, material orthotropy, and transverse load distributions are investigated. The classical plate theory, employing the Kirchhoff-Love hypotheses, is used in [68]. Simply supported and clamped-free boundary conditions on the unloaded edges are considered. Combined experimental and finite element techniques were suggested by Chai [69] to elucidate the postbuckling response of unilaterally constrained plates under monotonically increasing edge forces. Bifurcation phenomena for buckled partially and totally clamped plates resting on nonlinear elastic foundations are investigated in the work [38]. The influence of linear Winkler-type and shear nonlinear Pasternak-type foundation parameters is investigated.
A detail analysis of influence of Winkler-type foundation and other types foundations can be found in the works of [56,57].
Another well-known approach for constructing contact models is a penalty method. The variational principle with penalty is used for development of variational unilateral contact problem in the works of Ohtake et al. [32,33]. The variational model with penalty for buckling was also considered by Muradova and Stavroulakis [36]. In this case, Equation (6) becomes inequality Here, p is transverse reaction on the plate and p(w + b) = 0, w ≥ −b.
If we introduce the integral forms then Equations (2) and (12) are written as variational inequality, where the functions z, v ∈ W 2,2 (G) (W 2,2 is the Sobolev space). Applying the penalty method to Equation (13), one obtains where ε is a small positive number and π(w ε , z) is the penalty, defined as:

Delamination of the Buckled Plate
In this section, we discuss about mathematical contact models for buckled composite plates with delamination. The composite plates consist of several layers. After applying in-plane forces, different layers become separated (delaminated) and then buckle away from each other. Most commonly, a thin layer will buckle away from a thicker layer, typically of lower modulus of elasticity. When the bonding action between the layers is ignorable, this type of buckling problem can be considered in a state of unilateral contact.
The unilateral contact buckling behavior of delaminated rectangular plates in a composite member is analyzed, e.g., in the article of Ma et al. [70]. In the paper of Biggers and Pageau [71], tailored plates with buckling are investigated. Shen and Li [72] investigated postbuckling responses of shear deformable laminated plates resting on elastic foundation, subjected to in-plane compressive edge loads and a uniform temperature. The plates are constrained by an elastic foundation. The formulations are based on the large deflection plate theory, i.e., the von Kármán-type of kinematic non-linearity, and take into account the plate-foundation interaction. In the work [73], a delamination of buckled laminated plates is considered. Each laminate is modeled as an orthotropic Mindlin plate. For numerical solving, a combination of finite elements and asymptotic expansion is applied.
Uy and Bradford [74] studied numerically buckling of composite steel plates. Ovesy and Kharazi in [75] investigated the compressive bucking and postbuckling behavior of composite laminates with through-the-width delamination. The analytical method is based on the first-order shear deformation theory, and its formulation is developed on the basis of the Rayleigh-Ritz approximation technique.
The behavior of elliptical sublaminates, created by delaminations in composite plates, was investigated by Peck and Springer [76]. The plates are subjected to in-plane compressive, shear, and thermal loads. The developed model takes into account the stresses, strains, and displacements of the sublaminate as well as the loads applied to the plate.
The laminated composite plates are modeled mathematically using two higher-order shear deformation theories in conjunction with finite elements in the work of Sahoo et al. [77]. The final form of the governing equations of the bending and the free vibration responses are obtained using the variational method and the classical Hamilton principle, respectively.
In the paper of Yeh et al. [78], the Lagrangian formulation is proposed to analyze the bending behavior of the laminated plates and the local buckling phenomenon of the sublaminates in the delaminated region.
A mathematical model for analyzing delaminated plates is developed using three-dimensional (3D) degenerated elements in the work of Tiwari et al. [79]. The elements are introduced, applying the degenerated solid approach based on the Reissner-Mindlin assumptions.
In the article [80], a relevant survey on the various analytical models and numerical analysis for the free vibration of delaminated composites is provided. A basic understanding of the influence of the delamination on the natural frequencies and the mode shapes of composite laminates is considered.
Some delamination mathematical models are issued from the classical linear plate theories, namely from the Kirchhoff-Love and the Mindlin-Reissner principles, although many works related with buckling and delaminations of the plate are solved directly by FEMs or other numerical techniques (e.g., Obdrzalek and Vrbka [81]).
Below, we consider an example of unilateral contact buckling behavior of delaminated plates in a composite member. Buckling of two contacting plates is shown in Figure 4.  In [70] in order to describe a behaviour of two interacting plates, the linear deflection Kirchhoff-Love equation for each plate is used, i.e., where w ic and w i are the vertical displacements in the contact and non-contact regions, respectively, of plate i, D i is the rigidity of plate i, q r is the contact force between the two plates and σ xi is the normal stress in the x-direction. Delamination and unilateral contact in nonsmooth mechanics are investigated, e.g., by Stavroulakis and Panagiotopoulos in [28] and Demyanov et al [82]. A mathematical model for frictionless, adhesive and contact between a viscoelastic body and a rigid obstacle is studied, for example, in the article of Chau et al. [83].

Nonlinear dynamic equations
In comparison with the static case the dynamic initial-boundary value problem for the two coupled nonlinear dynamic equations has a unique solution according to the classical theory of nonlinear elasticity (Lasiecka [84] and Matkowsky and Putnick [85]).
The nonlinear plate dynamic equations are studied, for example, in the works of Shivamoggi [86], In [70], to describe a behavior of two interacting plates, the linear deflection Kirchhoff-Love equation for each plate is used, i.e., where w ic and w i are the vertical displacements in the contact and non-contact regions, respectively, of plate i, D i is the rigidity of plate i, q r is the contact force between the two plates, and σ xi is the normal stress in the x-direction. Delamination and unilateral contact in non-smooth mechanics were investigated, e.g., by Stavroulakis and Panagiotopoulos [28] and Demyanov et al. [82]. A mathematical model for frictionless, adhesive, and contact between a viscoelastic body and a rigid obstacle is studied, for example, in the article of Chau et al. [83].

Nonlinear Dynamic Equations
In comparison with the static case, the dynamic initial-boundary value problem for the two coupled nonlinear dynamic equations has a unique solution according to the classical theory of nonlinear elasticity (Lasiecka [84] and Matkowsky and Putnick [85]).
Several numerical models for treating the dynamic system have been developed, for example, Nath and Kumar [89] developed one that uses Chebyshev's series. Finite differences and C 1 h-version FEMs are employed for the dynamic nonlinear model by Gordnier and Fithen in [90]. Kirby and Yosibash [91] applied the pseudo-spectral method in space.
In comparison with the static case where we can only consider compressive forces, in the dynamic system, we can study both buckling and stretching phenomena together as well as separately. The nonlinear dynamic system of two coupled nonlinear hyperbolic type partial differential equations of vibrations of a rectangular isotropic elastic plate, taking into account viscosity, reads where w(t, x, y) is the deflection function of time and spatial variables (vertical displacement of the plate), ψ(t, x, y) is the Airy stress function of time and spatial variables, describing as before internal stresses appearing due to the deformation of the plate (e.g., [1,3]), and q(t, x, y) are the external disturbance forces. The domain of definition is Ω = (0, T] × G, where T is the final time and G is the shape of the plate as before. Further, in the system in Equations (14) and (15), ρ is the density of the material and c is the structural viscosity coefficient. The operator L t in Equation (14) is a function of time and characterizes an external (compressive or tensile) constant force applied to the edges of the plate. It is defined according to the above-mentioned assumptions (see [48]), Here, λ = λ(t) and λ = λ (t) are the factors of moving compressive and tensile loads uniformly distributed on the edges of the plate. If q(t, x, y) = 0, λ = 0, and λ = 0, then we have free vibrations.
Equations (14) and (15) describe vibrations of the plate, subjected to compressive and tensile (stretching) non-constant loading, applied simultaneously or separately. The exciting forces can change with time and they are uniformly applied on the edges of the plate. If the Monge Ampère part does not present, then we have the dynamic linear plate model.
If the structure is also unilaterally supported by the upper and/or lower elastic foundations with different stiffnesses, then, as discussed in Section 3, the transversal loading q = −p(w, ...), i.e., the dynamic buckling problem, involves contact phenomenon. The cases when the nonlinear foundations are modeled in terms of the nonlinear elastic Winkler-type and shear Pasternak-type are considered in the paper of Muradova and Stavroulakis [39].
In the work of Kanto and Yagawa [92], the penalty function method is applied to incorporate the contact conditions in the equation of motion and a trial and error method is used. The technique is applied to a dynamic buckling problem involving contact phenomenon.

Conclusions
Mathematical models with buckling and contact phenomena are discussed in this paper. The effects of buckling, which occurs when in-plane forces are applied on the edges of the plate, is shown on examples. Different types of boundary conditions for the plate are discussed.
The contact phenomenon is caused by the presence of elastic foundations or interaction of parts of a delaminated composite plate. The contemporary models with buckling, including subgrade model for linear and nonlinear foundations, penalty method, used in studies with buckling, namely subgrade model, and delamination plates are described. The semi-analytical and numerical techniques for resolving these problems are also briefly discussed. The references relevant to the considered topics are given. The advantages and disadvantages of the previous research are shown.
The buckling and contact problems have a long history. The basic principles of understanding of buckling and contact effects arising in elastic plates and the methods of treating the problems related with them are discussed in the previous sections. Based on the previous research results, one can conclude that the mathematical models, combining buckling and contact phenomena together are studied, in the majority of the publications, either with the help of finite element methods or by using analytical or semi-analytical approaches [93,94].
The finite element analysis is suitable for the mechanical problems for plates with complicated geometry and contact nonlinearity. It can provide meaningful and accurate results regardless material properties, boundary conditions, and loading. FEMs also work perfectly for bifurcation analysis (e.g., [42,44,45,49]) and contact and delamination problems (e.g., [62,66,69,73,77]). Together with variational Galerkin method, they can provide effective analysis for buckling nonlinear plate problems [42]. A mixed Galerkin-perturbation technique can be used to determine thermal buckling loads and post-buckling equilibrium paths. Effective mathematical study of different types of buckling such as plastic, dynamic, and external is supported by FEMs [94].
Although FEMs are most widely used and dependable in buckling models, contact problems, and buckling analysis of laminated plates plate problems, there are some disadvantages in applications of them. In particular, FEMs are sensitive to mesh density, the rate of convergence in simulation, element type, computational time, and shear locking. Shear locking is an error that occurs in finite element analysis due to the linear nature of quadrilateral elements.
FEMs are flexible for the geometry of the structure, but analytical and semi-analytical methods provide higher accuracy. FEMs require large amount of computational resources for reaching high accuracy, which may be more suitable for detailed analysis rather than overall demonstration of structural performance. In contrast, analytical or semi-analytical approaches provide a high accuracy of computations and can show overall behavior of plates (e.g., [88,95]).
Most of the analytical buckling approaches are based on the energy principles ( [82,93]). Approximation models for buckling problems usually employ energy principle or semi-analytical approaches to determine the critical buckling load of structures, of which the Rayleigh-Ritz method, Kantorovich method, finite strip method, etc. are commonly used (e.g., [75,93]).
Note that delamination usually involves contact leading to nonlinear post-buckling problems, and analytical models might not give consideration of all the factors simultaneously. Numerical simulations and experiments are needed to provide further investigation on the phenomenon.
The disadvantages of the analytical methods are their restricted applicability to relatively simple shapes and boundary conditions, as well as to buckling and contact of plates with finite sizes. Some assumptions and large simplifications should be done in material model, geometry, and boundary conditions. These obstacles can be decreased by applications of semi-analytical methods for some cases. For example, for rectangular elastic plates, which is widely used in civil engineering, the semi-analytical methods based on variational principles (spectral or pseudo-spectral methods) work best and can provide a high accuracy of computations ( [38,47,89,91]).
Summarizing, one can say that the choice of the method for treating the coupled buckling and contact problems in elastic plates depends on many factors, such as material properties, boundary conditions, required accuracy of computations, geometry of the structure, etc. The advantages and disadvantages of the techniques should be taken into account.
Although sufficient research work has been done for the investigation buckling and contact phenomena together for thin plates, there are still open questions in creation of new effective models or modification of the existent ones and then applications of them to complex nonlinear problems in science and engineering.
Investigations of the nonlinear systems including buckling and contact effects are keys for solving various problems in biomechanics, cell wall mechanics, and deformations of thin membranes, where exploration of interactions between buckled plate and elastic foundation or buckled parts of composite structures is important and challenging. Therefore, since most nonlinear models cannot be solved analytically, it is necessary also to develop effective semi-analytical and/or numerical techniques allowing to overcome barriers in simulation of complex systems.