## Abstract

We present a comprehensive study on the postbuckling response of nonlocal structures performed by means of a frame-invariant fractional-order continuum theory to model the long-range (nonlocal) interactions. The use of fractional calculus facilitates an energy-based approach to nonlocal elasticity that plays a fundamental role in the present study. The underlying fractional framework enables mathematically, physically, and thermodynamically consistent integral-type constitutive models that, in contrast to the existing integer-order differential approaches, allow the nonlinear buckling and postbifurcation analyses of nonlocal structures. Furthermore, we present the first application of the Koiter’s asymptotic method to investigate postbifurcation branches of nonlocal structures. Finally, the theoretical framework is applied to study the postbuckling behavior of slender nonlocal plates. Both qualitative and quantitative analyses of the influence that long-range interactions bear on postbuckling response are undertaken. Numerical studies are carried out using a 2D fractional-order finite element method (f-FEM) modified to include a combination of the Newton–Raphson and a path-following arc-length iterative methods to solve the system of nonlinear algebraic equations that govern the equilibrium beyond the critical points. The present framework provides a general foundation to investigate the postbuckling response of potentially any type of nonlocal structure.

## 1 Introduction

The study of buckling and postbuckling responses of elastic structures have obvious applications in the domain of structural analysis and design. Seminal work in the study of postbuckling response of continuum structures was conducted by Koiter [1]. Later, Budiansky [2] extended this study to the postbuckling response of elastic structures. These works generated broad interest in the academic and research communities involved in the structural design [3]. By employing models and methodologies derived, either directly or indirectly, from these seminal works, the postbuckling response of several structures have been successfully studied [4–8]. These studies are based on the classical (local) constitutive theories that involve a point-wise correspondence between the displacement and the stress field generated at the point. Although the localized relations of constitutive theories are still fundamental for most studies on the behavior of solids, the existence of nonlocal interactions between distant points in a solid have been established through experimental investigations. These observations point toward a nonclassical behavior that is not accounted for in classical local elasticity theories. Some common examples of solids exhibiting nonlocal behavior include (microscale and nanoscale) structures where size effects gain prominence [9,10]. In addition, recent experimental observations point to a more pervasive presence of nonlocal interactions even in macroscale applications such as sandwich structures, stiffened panels, and functionally graded materials [11,12].

Owing to the importance of the applications mentioned earlier, several theories have been proposed to model the nonlocal effects in elastic solids. Among the numerous theories existing in the literature, the strain-based integral constitutive relations for the nonlocal solids proposed by Kröner [13] and Eringen [14] are likely the most prominent. In these approaches, nonlocal interactions are introduced in the constitutive law by means of an integral stress-strain relation defined over an assigned domain of influence. This integral effectively captures long-range interactions by means of a convolution of the strain tensor weighted by a proper attenuation function. While this approach has certain key advantages, primarily in providing an intuitive representation of the long-range interactions via convolution integrals, it suffers some key mathematical and physical inconsistencies. The salient deficiency of integral models is that it belongs to an ill-posed class of integral equations that do not admit unique solutions [15–17]. To circumvent the implicit formulations associated with the integral approaches, an equivalent gradient-based model of nonlocal elasticity was developed [18]. The modeling advantages introduced by these gradient-based approaches sparked numerous investigations focused on assessing the effect of long-range interactions on the response of structural elements [19–21]. The equivalence between the gradient and integral models is based on certain restrictions like the nonlocal influence domain being unbounded and on certain choices of the attenuation kernel. These conditions place physical restrictions that do not allow the use of gradient models for a general study on nonlocal elasticity. This is because violation of these restrictions would present paradoxical observations, such as the vanishing nonlocal effects, encountered for certain choices of loading and boundary conditions [15,17]. In the context of stability analysis, the nonlocal effects are realized by a modification of the system stiffness terms. When employing integral-type nonlocal constitutive laws, the material stiffness is reduced due to nonlocal interactions [22,23], whereas differential models cause an increase in the geometric stiffness [18,24]. In either case, both models predict a lower critical load in nonlocal structures when compared with their classical local counterpart [23]. Note that most of the literature addressing the effects of nonlocal mechanisms on the critical load has focused on linear buckling. Studies on the critical load in the case of nonlinear buckling and postbuckling response are almost entirely based on differential models. Among the existing studies, nonlocal effects are restricted to carbon nanotubes (CNTs) and graphene. Both classical and modified differential models of nonlocal elasticity were employed to compare the effects of long-range interactions on both the hardening and softening response of low-dimensional structures [25–27]. Additional studies include postbuckling of nonlocal beam and cylinder models for CNTs [28,29] and plate and shell models for graphene sheets [30,31]. Nonlocal effects have also been considered in the study of the postbuckling response of bio-based protein tubules [32] and in the pull-in instability of the micro- and nano-electromechanical devices [33,34]. To the best of author’s knowledge, all the reported studies are based on the differential models of nonlocal elasticity and conclusively predict a lower value of critical load across boundary conditions followed by a reduction in the postbuckling resistance to deformation.

The lack of similar studies based on integral models of nonlocal elasticity may be attributed to the relative complexity of the integral-type constitutive relations, and the unavailability of a consistent and positive semi-definite definition of the deformation energy density functional within this framework. In fact, it has been argued that such a functional cannot be defined for the nonlocal solids [16]. Since postbuckling analysis is based on energy approaches, this latter point presents a significant technical gap that requires to be addressed. Moreover, it must be noted that the integral and differential-type constitutive laws discussed earlier are thermodynamically inconsistent [35]. This is because of the ad hoc constitutive laws that do not satisfy the thermodynamic balance in a rigorous manner. More specifically, when nonlocal solids are modeled via integral or differential approaches, the first and second law of thermodynamics are only satisfied in a weak manner (over the nonlocal domain) as opposed to the physically consistent localized (strong) form [35]. Thus, a consistent energy formulation may not be developed employing these models. This consideration extends also to the alternative integral-type two-phase definition of constitutive law that includes both local and nonlocal terms in the stress–strain relations. While the two-phase definition admits unique solutions by virtue of the well posedness of Fredholm’s second integral, it does not satisfy the second law of thermodynamics (Clausius–Duhem inequality) in a rigorous manner. This highlights the limitations of the existing constitutive models for nonlocal elasticity and their deficiencies in an analysis over postbuckling studies [36].

Fractional calculus has demonstrated enormous potential for addressing areas where integer-order models fall short, such as the case of modeling of complex phenomena. Limiting the discussion to the current context of structural models, recent developments include fractional-order constitutive theories of nonlocal elasticity that can potentially address the aforementioned deficiencies [37–41]. Other prominent examples of successful applications of fractional-order theories are the constitutive modeling of memory effects in viscoelastic materials [42,43] and nonlocal effects in a multiscale analysis [38,41,44]. In the case of fractional-order constitutive theories, the differ-integral nature of fractional-order derivatives serves as a suitable alternative to the integro-differential constitutive relations encountered in nonlocal elasticity. Indeed, the fractional-order derivatives appearing in the kinematic strain–displacement relations allow accounting for long-range interactions within the constitutive model of a nonlocal elastic solid. Successful applications of this theory to the development of mathematical and physically consistent models for softening [41] and stiffening [44] type of nonlocal interactions have also been proposed. The frame-invariant fractional-order continuum theory based on a physically consistent definition for position-dependent nonlocal length scales is particularly effective in allowing a rigorous application of the thermodynamic principles [36]. This observation is crucial in developing consistent constitutive laws for nonlocal elasticity because it allows defining a positive semi-definite functional corresponding to the deformation energy density of a nonlocal solid. It has been proved that the positive semi-definite form of the deformation energy density functional leads to self-adjoint linear operators in the governing equations [45,46]. By employing weak-form expressions and variational methodologies, finite element techniques suitable for the numerical solution of boundary value problems based on fractional-order equations have also been developed [45].

Following energy-based methods afforded thanks to the fractional-order model, softening effects due to nonlocal interactions have been accounted for both the linear and nonlinear response of beams and plates [45–48]. The large deformation analysis of the nonlocal solid requiring geometrically nonlinear strain–displacement relations is particularly insightful within the context of stability of fractional-order continuum theory. Unlike the integral and differential models of nonlocal elasticity that modify either the material or the geometric stiffness term alone, physically the nonlocal interactions are expected to be realized on both of these stiffness terms. This is achieved by incorporating the nonlocal interactions at the level of kinematic relations within the constitutive relations for a nonlocal solid. Thus, the effect of the nonlocal interactions over the critical load is decided by the net result of modifications to the material and geometric stiffness terms by the fractional-order derivatives. This observation has been drawn following rigorous parametric studies isolating the nonlocal interactions in each term to identify the net result over the critical load corresponding to linear buckling [49]. It is imperative that a study of the postbuckling response of nonlocal structures, which also depends on both of these stiffness terms, is conducted employing the mathematically, physically, and thermodynamically consistent nonlocal model of the fractional-order continuum theory.

With this objective in mind, we study the critical load corresponding to the nonlinear buckling and the subsequent postbuckling response of nonlocal structures. For this purpose, we build upon the existing framework for the stability analysis of fractional-order structures [49]. As mentioned earlier, this study is possible due to the positive semi-definite deformation energy density functional made available by the fractional-order formulation of nonlocal elasticity [36,45]. The objective of the current work is twofold. First, the nonlinear buckling of fractional-order structures is studied by identifying the singular points of the equilibrium curves. We develop the nonlinear equations required to obtain the critical load and the associated buckling mode following variational principles. Subsequently, the postbifurcation analysis in the immediate vicinity of the singular points is conducted. Second, the developed framework is employed to perform a numerical study (via fractional-order finite element method (f-FEM)) and to assess the effect of nonlocal interactions on the postbuckling response of a slender plate.

This article is organized as follows. First, it provides a brief summary of the fractional-order continuum theory. This is followed by a summary of the energy-based framework employed for the identification of critical points of nonlocal structures corresponding to linear buckling. Then, it presents the derivation of the critical loads corresponding to singular points of nonlinear equilibrium. Thereafter, by means of extending the Koiter’s asymptotic method, the postbuckling response in the immediate vicinity of the critical points is studied. Finally, an analysis of the postbuckling behavior of fractional-order Kirchhoff plates is performed. We provide brief details regarding the 2D f-FEM model developed for this numerical study over the postbuckling response of nonlocal structures.

## 2 Postbuckling Response of Fractional-Order Nonlocal Structures

In this section, we develop the models required to study the postcritical response of fractional-order nonlocal structures. To provide the necessary mathematical background and notation, we briefly review the basic constitutive relations for the fractional-order continuum theory [41,47] and the stability of equilibrium points for a nonlocal solid evaluated following this theory [49]. We use the boldface notation to describe tensors (including vectors), while subscript indicates the measure number of the component evaluated in Cartesian coordinates.

### 2.1 Constitutive Model for Fractional-Order Structures.

**U**(

**X**) is the Lagrangian displacement field. The term $\u2207\alpha U$ denotes a Riesz–Caputo (RC) fractional derivative of order

*α*∈ (0, 1] of the displacement field with respect to the spatial coordinates $X\u2208R3$. This space-fractional derivative is defined as follows (no summation over repeated indices) [45]:

*U*

_{i}(

**X**) with respect to the coordinate

*X*

_{j}. The terminals of the Riesz–Caputo fractional-order derivative defined at

*X*

_{j}are $XAj$ and $XBj$. These terminals along the

*X*

_{j}–direction are obtained using position-dependent length scale parameters $LAj$ and $LBj$ as $XAj=Xj\u2212LAj$ and $XBj=Xj+LBj$. Thus, the domain $(XAj,XBj)$ defines the horizon of nonlocal interactions along

*j*-direction at the point

**X**. Finally, Γ(·) is the Gamma function, and the term 1/2Γ(2 −

*α*) along with the position-dependent length scales allow the second-order fractional-order strain tensor to be frame invariant [41]. The above given RC definition for the fractional derivative may be recast using the definitions of the Caputo derivatives as follows:

*α*-order power-law kernel that serves as the attenuation function for nonlocal interactions between points

**X**and $\xi $. The differ-integral nature of the fractional-derivative evident from above result introduces the long-range nonlocal interactions into the constitutive models at the level of the kinematics.

**U**(

**X**), respectively. Note that, following the assumptions for the von-Kármán displacement relations, the quadratic component includes the fractional-order derivatives of only the transverse displacement field

*U*

_{3}(

**X**) with respect to in-plane coordinates. The transverse strains (normal and shear) are simply the linearized forms of complete nonlinear strain given in Eq. (1), which is $e~ij(U)$ (see Eq. (4)). For the complete expressions of the nonlinear fractional-order von-Kármán strain–displacement relations, the reader is referred to Ref. [47]. Following this review over the fractional-order strain–displacement relations for a nonlocal solid, we will discuss the stress–strain laws corresponding to the current formulation.

*λ*and

*μ*. Finally, the localized material constitutive relations corresponding to the nonlocal response of a fractional-order solid may be written as follows [36]:

### 2.2 Critical Points of Stability for Fractional-Order Structures.

**u**(Λ); Λ] of this nonlocal structure is expressed as follows:

**f**(Λ) applied on the boundary ∂Ω

^{σ}. The first variation of the aforementioned expression for the potential energy

*δ*Π[

**u**] = Π,

_{u}

*δ*

**u**provides governing equations of equilibrium for the fractional-order structure [45,49]. The configuration given by the coordinates

**u**

_{e}(Λ

_{e}) is an equilibrium state for the nonlocal structure if it satisfies the equation

*δ*Π[

**u**

_{e}] = 0 given by

*δ*

**u**. In the derivation of the aforementioned equation, we make use of the expression for the deformation energy density of a nonlocal solid given in Eq. (5) and the von-Kármán strain–displacement relations in Eq. (4). Note that the comma notation in the subscript used earlier indicates a differentiation, but not necessarily to an integer order. For further details regarding the variational principles and the transversality conditions for fractional-order systems, the reader can refer to Ref. [50].

**u**

_{e}(Λ

_{e}) is stable if the second-variation

*δ*

^{2}Π[

**u**

_{e}] = Π,

_{uu}(

**u**

_{e})

*δ*

**u**

_{1}

*δ*

**u**

_{2}> 0, where

*δ*

**u**

_{1}and

*δ*

**u**

_{2}are independent variations [49]. Before we proceed further, some assumptions should be mentioned. Note that these assumptions are common to similar studies on bifurcation of elastic solids using classical elasticity [2]. We assume a small predeformation at the critical point. By this, we consider an equilibrium curve

*C*

_{0}for Λ < Λ

_{c}given by

**u**=

**u**

^{0}(Λ), such that

**u**

^{0}≪

**v**. Here,

**v**is the postbifurcation response of the nonlocal structure. Furthermore, we assume a proportional loading force (

**f**(Λ) = Λ

**f**

^{0},

**f**

^{0}being a representative force vector) that results in the stress generated to be given as $\sigma ~(\Lambda )=\Lambda \sigma ~0$. We identify the critical configuration of stability

**u**

_{c}(Λ

_{c}) as the solution to the variational equation defined over

*δ*

^{2}Π[

**u**

_{c}] as follows:

**U**

_{b}=

*δ*

**u**

_{1}is the critical mode corresponding to the bifurcation at

**u**

_{c}and Λ

_{c}is the critical value of the load control parameter. The configuration (

**u**

_{c}, Λ

_{c}) that is a solution to the aforementioned nonlinear equation provides the singular point of the nonlinear elastic response of the nonlocal solid. For ease of notation, we express

*δ*

^{2}Π[

**u**

_{c}] evaluated at

**u**

_{c}(Λ

_{c}) as

*δ*

^{2}Π

^{c}. The aforementioned equation may be further simplified through linearization to derive the following expression for the critical points corresponding to

*linear*buckling of fractional order solids [49]:

We point to a simultaneous modification of the standard *and* geometric stiffness of the nonlocal structure (i.e., the numerator and denominator in the above expression, respectively) caused by the long-range interactions. This is because of the fractional-order derivatives being present in the expressions of both stiffness terms in the aforementioned result for critical load of nonlocal structure. Therefore, the effect of the nonlocal interactions on the critical load $\Lambda cl$ is not trivial and may either increase or decrease due to long-range effects. This is in contrast with the observations following from existing integer-order nonlocal theories that predict a consistent decrease in the critical load. We clarify this statement by a comparison of the fractional-order model with the standard integer-order Eringen’s integral and differential theories of nonlocal elasticity. Following the integral theory [14], the effect of the nonlocal interactions is captured only within the material stiffness terms. This causes a reduction of the material stiffness with an increasing degree of nonlocality. However, the geometric stiffness evaluated using this theory is identical to the corresponding result following from classical (local) elasticity [23]. Conversely, the Eringen’s differential theory increases the geometric stiffness terms as a consequence of the nonlocal interactions, while leaving the material stiffness unchanged from classical (local) elasticity predictions [24]. Please refer to Ref. [49] for a detailed comparison of the critical buckling load evaluated for a nonlocal solid employing the fractional-order and integer-order nonlocal continuum theories. While these observations are drawn for the linear buckling of fractional-order structures, they are relevant for nonlinear buckling studied using the nonlinear eigenvalue problem given in Eq. (9).

### 2.3 Postbifurcation Analysis of Fractional-Order Structures.

Following the previous study on the identification of singular points of the fundamental solution of the governing equations, our next objective is to analyze the nature of the singularity and develop the models required for the subsequent analysis of the nonlinear response of fractional-order solids beyond these critical points. In the domain of classical elasticity, this study of the postbuckling response is carried out following the energy methodology [1]. Here, we make use of the positive semi-definite deformation energy density (see Eq. (5)) achievable via the fractional-order derivatives for nonlocal solids. We proceed with a straightforward extension of the Koiter’s asymptotic method, originally defined for integer-order models of classical elasticity, to obtain the postbuckling response of fractional-order nonlocal structures.

**u**, Λ) adjacent to the critical point (

**u**

_{c}, Λ

_{c}) in terms of the parameter

*χ*[2]:

**v**is the postbifurcation response expressed as the difference between the fundamental curve

**u**

_{c}∈

*C*

_{0}(

**u**

^{0}) and the bifurcated curve

*C*

_{1}(

**u**). This difference is expressed in terms of the critical mode

**U**

_{b}obtained in Eq. (9) and correction factor

**W**defined such that

**W**⊥

**U**

_{b}. The quantities (

**U**

_{b}, Λ

_{I}) and (

**w**

_{2}, Λ

_{II}) provide the slope and curvature of the load–displacement equilibrium curves at the point

*χ*= 0 [51]. As in the case of classical elasticity, in the immediate vicinity of the critical point, a linear approximation would suffice. In case the coefficient of the linear approximation (Λ

_{I}) is zero, the higher order terms will be considered. Assuming the potential energy functional to be continuous in the vicinity of the above identified critical point, the coefficients of the aforementioned expansion can be obtained analogous to integer-order models [1,2]. It must be mentioned that by making use of the Taylor series expansion about the critical point (

**u**

_{c}, Λ

_{c}), the error in approximation approaches zero for Λ → Λ

_{c}. Employing the asymptotic approach to the potential energy expansion for nonlocal structures, our objective is twofold: (1) establish the bifurcation of the primary equilibrium curve at the critical point and (2) discuss the stability of the bifurcated equilibrium curves beyond the equilibrium. Here, for the sake of brevity, we shall discuss only the salient conclusions regarding the influence of the nonlocal interactions drawn from this mathematical exercise.

*C*

_{0}(when the control parameter Λ increases from 0), the second variation of the potential energy

*δ*

^{2}Π is positive. This corresponds to the stability of equilibrium positions along this curve. However, for increasing values of Λ,

*δ*

^{2}Π reduces until finally at the critical point Λ

_{c}we have

*δ*

^{2}Π

^{c}= 0. Therefore, assuming the independent variation

*δ*

**u**

_{2}=

**U**

_{b}(critical mode shape) in Eq. (9), we obtain [2]:

*C*

_{0}undergoes an angular bifurcation at the critical point (

**u**

_{c}, Λ

_{c}) in a nonlocal solid. This condition is analogous to the one known in classical elasticity where angular bifurcation is noted for the nonlinear equilibrium curves at the critical (singular) point. From the aforementioned result, we may also conclude that the nonlocal interactions captured by the fractional-order derivatives do not influence the nature of the bifurcation. This is because the quadratic term $q~(v,v)>0$ for nonzero

**v**. This conclusion is necessary to ascertain the nature of the instability in nonlocal solids at the critical point. While studies on postcritical response of nonlocal structures are available in the literature [28,30,31], these works based on the integer-order models of nonlocal elasticity assume angular bifurcation without any formal proof. Following this discussion, we now prove the nature of the instability in the presence of nonlocal interactions. This result is possible due to the positive semi-definite quadratic term $q~(v,v)$ available from the fractional-order kinematic relations for nonlocal solids. This observation gains further significance in the analysis of postbuckling response of complex nonlocal structures with a priori unknown nature of the instabilities. In accordance with this result, the conditions required for instability in nonlocal solids simply follow from classical elasticity [1]. First, the nonzero quadratic term following from the von-Kármán strain–displacement relations (geometric nonlinearity) is essential to ensure a nontrivial result in the aforementioned expression. Second, a compressive force

**f**

^{0}causing $\sigma ~0<0$ is also required to induce an instability within the nonlocal structure.

*C*

_{0}and the bifurcated curve

*C*

_{1}beyond the critical point. To draw a qualitative conclusion, we look at the sign of the difference in potential energy between the two configurations on the

*C*

_{0}and

*C*

_{1}curves ΔΠ = Π(

**u**, Λ) − Π(

**u**

^{0}, Λ) for a given value of the load parameter Λ. This result for the nonlocal structure can be derived following the asymptotic relations given in Eq. (11) and the transversality rate defined in Eq. (12) to be as follows:

*a*

_{1}and

*a*

_{2}are positive constants. For $T<0$, we have ΔΠ < 0 (see Eq. (12)). Thus, equilibrium configurations (

**u**, Λ) along the bifurcated curve

*C*

_{1}admit lower potential energy when compared with the fundamental branch for Λ > Λ

_{c}. We conclude that the bifurcated branch at the critical point Λ

_{c}presents stable equilibrium solutions for Λ > Λ

_{c}.

Increasing the degree of nonlocality would result in lower values for fractional derivatives. Thus, the numerical values of $\sigma ~0$ and $q~$ in Eq. (10) will be reduced for the fractional-order structure. This would result in lower value of $T$ for fractional-order solids with the increasing degree of nonlocality. Therefore, we infer from Eq. (13) that nonlocal effects will reduce the difference in energy barrier between the fundamental and the bifurcated paths at the singular points.

## 3 Postbuckling Analysis of Fractional-Order Kirchhoff Plates

In this section, we apply the aforementioned formulation to study the postbuckling response of a nonlocal Kirchhoff plate. We begin with the derivation of the governing equations of equilibrium for a fractional-order Kirchhoff plate following variational principles. We also develop the numerical model for these nonlinear fractional-order (integro-differential) governing equations employing the path-following method based on arc-length constraints.

### 3.1 Basic Formulation for Fractional-Order Kirchhoff Plates.

*a*×

*b*, and the thickness is denoted by

*h*. The schematic of the plate in Fig. 1(a) indicates the Cartesian coordinates chosen for the current study. As shown in this figure,

*x*

_{3}= 0 is the mid-plane and

*x*

_{3}= ±

*h*/2 denote the top and bottom surface of the flat plate. As indicated in the schematic, the vertical free surfaces of the plate coincide with

*x*

_{1}= 0,

*a*and

*x*

_{2}= 0,

*b*. Following the slender assumptions of the Kirchhoff plate theory, the height of the plate is chosen such that

*h*<

*a*/50. Under this assumption, the components of the displacement field distribution

**u**(

**X**) at a point

**X**(

*x*

_{1},

*x*

_{2},

*x*

_{3}) ∈ Ω according to the Kirchhoff plate theory are given by

*u*

_{0}(

**X**

_{0}),

*v*

_{0}(

**X**

_{0}), and

*w*

_{0}(

**X**

_{0}) are the generalized displacement coordinates corresponding to the displacements at the point

**X**

_{0}(

*x*

_{1},

*x*

_{2}) on the reference plane

*x*

_{3}= 0. The fractional-order geometrically nonlinear strains evaluated using the von-Kármán strain–displacement relations in Eq. (4) are now given as follows:

*a*)

*b*)

*c*)

*d*)

*α*∈ (0, 1) along the

*x*

_{1}– and

*x*

_{2}–directions are denoted by $Dx1\alpha \u2261x1\u2212LA1Dx1+LB1\alpha $ and $Dx2\alpha \u2261x2\u2212LA2Dx2+LB2\alpha $. The length scales $LAi$ and $LBi$ (

*i*= 1, 2) in the aforementioned expressions indicate the nonlocal horizon of influence in

*x*

_{i}-direction at the point of interest

**X**

_{0}(

*x*

_{1},

*x*

_{2}). The RC derivatives model long-range interactions at the point

**X**

_{0}(

*x*

_{1},

*x*

_{2}) over the domains $(x1\u2212LA1,x1+LB1)$ and $(x2\u2212LA2,x2+LB2)$.

*a*)

*b*)

*E*is the Young’s modulus and

*ν*is the Poisson’s ratio for the isotropic solid.

Following the Kirchhoff plate theory, we note that the in-plane stiffness of the plate is much larger than the bending stiffness [52]. This results in transverse displacements, perpendicular to the mid-plane, being considerably larger when compared with the in-plane displacements. Furthermore, transverse shear rigidity of the slender plate is at least an order of magnitude higher than the bending rigidity. Therefore, the contribution of nonzero transverse shear strains to the deformation energy density of the fractional-order Kirchhoff plate may be neglected [45,48].

**F**

_{1}= ±

*F*

_{1}

**e**

_{1}and

**F**

_{2}= ±

*F*

_{2}

**e**

_{2}are applied at the free-edges

*x*

_{1}= 0,

*a*and

*x*

_{2}= 0,

*b*as indicated in the schematic shown in Fig. 1(b). In keeping with our discussion in Sec. 2, these loads are expressed in terms of a control parameter Λ and unit forces $Fi0=1ei$ (

*i*= 1, 2). In addition to these in-plane loads, we also apply an additional load

**F**

_{3}=

*F*

_{3}

**e**

_{3}transverse to the mid-plane of the plate. This load is applied to activate the stable branch of the postbuckling nonlinear response path at the critical point [53,54]. The details of the eccentric load will be provided in the following. Note that the eccentric load (

**F**

_{3}) serves a role similar to that of the geometric imperfections used in post-buckling analysis [55,56]. Furthermore, the distribution of this transverse load is chosen to yield the mode shapes corresponding to linear buckling modes [1,2].

*δ*Π = 0.

### 3.2 Numerical Modeling.

*a*] × [0,

*b*] of the plate is divided into

*N*

_{e}uniform finite elements denoted by $\Omega ei$, such that $\u222ai=1Ne\Omega ei=\Omega $ and $\Omega ej\u2229\Omega ek=\u2205\u2200j\u2260k$. In the current study, the unknown generalized displacement coordinates are given by

*U*} may be expressed in terms of the corresponding nodal coordinates as follows:

*a*)

*b*)

*δ*Π = 0, the algebraic equations of equilibrium are given as follows:

*K*

_{S}], the global assembly of the nodal displacement coordinates {Δ}, and the force matrices {

*F*} are available elsewhere [48]. Note that the stiffness matrix in the aforementioned expression is nonlinear and depends on the configuration of the deformed system.

*F*

^{0}} is the representative force vector and Λ is an additional variable that serves as the load control parameter. The nonlinear algebraic equations given in Eq. (21) can be recast as follows:

*R*} is the vector of residual nodal forces for the fractional-order structure as defined in Refs. [47,48]. Based on the Newton–Raphson method, the increment (

*δ*{Δ},

*δ*Λ) in equilibrium position is governed by the following equation:

*a*)

*K*

_{T}] is the

*tangential*stiffness matrix. An additional constraint over the increments

*δ*{Δ} and

*δ*Λ:

*b*)

*δ*Λ = 0 and

*δ*{Δ} = 0 correspond to load–control and displacement–control increments, respectively.

## 4 Results and Discussion

In this section, we report the numerical results for the models developed in Sec. 3 and provide a quantitative assessment of the effect of nonlocal interactions over the postbuckling response of Kirchhoff plates. More specifically, the effect of the fractional order and of the nonlocal length scales on the nonlinear load–displacement curves will be specifically addressed. In all simulations, we choose an isotropic solid with Young’s modulus *E* = 1.09 MPa and the Poisson’s ratio *ν* = 1/3 [48]. The Young’s modulus chosen here applies to a class of soft materials that exhibit large deformations, which is essential for the intended study. However, note that the framework developed here is not specific to this particular choice of material properties.

The constitutive parameters corresponding to fractional-order continuum theory, which are the fractional-order *α* and the nonlocal length scale *l*_{f}, are not fixed a priori. They will be varied to identify the effects of these individual parameters on the postbuckling response. Therefore, their value will be provided wherever required. The nonlocal length scales for the isotropic solid are chosen such that for points sufficiently away from the boundaries we have $LA1=LB1=LA2=LB2=lf$. This symmetry is broken for points close to the boundaries of the plate. As illustrated in Fig. 2, for points whose distance from the boundary (in a given direction) is smaller than *l*_{f}, we truncate the nonlocal length scale at the physical boundaries of the plate. Concerning the numerical simulations carried out here below, we consider a square plate with in-plane dimensions *a* = *b* = 1 m. The thickness of the plate is chosen to be *h* = *a*/120 consistently with the assumptions for the Kirchhoff plate displacement theory.

*F*

_{1}applied at the edges

*x*

_{1}= 0,

*a*. As discussed in Sec. 3, we also apply a uniformly distributed transverse load to activate the stable branch of the postbuckling path at the critical point. Furthermore, the numerical results for two choices of boundary conditions will be analyzed. In the first case, all the edges of the plate are simply supported. The unloaded edges in this case are free to move laterally along the plane of the plate. This choice referred to as the SSSS-01 is given by [56,59]:

*a*)

*b*)

*v*

_{0}≠ 0) when subject to compressive forces on the two edges

*x*

_{2}= 0,

*b*. Also, in both cases, the plates are loaded on a movable edge (

*u*

_{0}, ∂

*w*

_{0}/∂

*x*

_{1}≠ 0). We further reiterate that the framework is general, but we restricted the numerical investigations to the above standard choices of boundary conditions [56,59]. All the numerical results reported here are nondimensionalized as follows:

### 4.1 Validation.

Before proceeding with the analysis of the postbuckling behavior, we inspect the efficacy of the numerical model that was developed. Specifically, we establish the convergence of the 2D f-FEM and then validate the numerical model by comparison with existing results available in the literature.

To establish the convergence of the 2D f-FEM, we compare the postbuckling responses obtained from successive refinements of the finite element mesh. As seen in other studies on nonlocal elasticity, the numerical approximation of the system matrices includes an additional integration [60]. In the case of the fractional-order strain–displacement relations employed here, this is the convolution integral present in the definition of the nonlocal strain (see Eq. (3)). The accuracy of the numerical approximation of this convolution integral is controlled by the “dynamic rate of convergence,” which is defined as $N=lf/le$. For the 2D f-FEM with convolution integrals defined in both the *x*_{1} and *x*_{2}-directions, we have two parameters in $Ni$ (*i* = 1, 2) that control the element discretization in the *x*_{i}-direction. The nonlinear load–displacement curves corresponding to different choices of this parameter are compared for the SSSS-01 plate in Fig. 4(a). From these results, we conclude that an excellent convergence of the postbuckling curves to within 1% difference is achieved for successive refinements beyond $N1\xd7N2=12\xd712$. A similar conclusion may be drawn for the convergence of the postbuckling curves for CSCS-01 plate from Fig. 4(b). On the basis of previous convergence results, we use $N1\xd7N2=12\xd712$ for all the numerical simulations presented in the following section.

Next we validate the postbuckling behavior obtained from the current model when using *α* = 1.0 (i.e., the local case) with the corresponding results in the literature on classical isotropic plates. The efficacy of the 2D f-FEM as a numerical tool for fractional-order plates was established in our previous studies [46,48], so we validate the numerical method and the efficacy of the arc-length constraints directly on a study of postbuckling response. Note that, for the choice of *α* = 1.0, the fractional-order continuum theory identically reduces to classical elasticity theory [41,45]. The postbuckling curves for the isotropic square plate with *α* = 1.0 and subject to SSSS-01 boundary conditions are shown in Fig. 5(a). In this figure, the nonlinear load–displacement curves are compared for different amplitudes of the eccentric load $F\xaf3$. We validate these curves with the corresponding result available from the classical theory [56].

A similar exercise is carried on for the validation of the postbuckling curves for the CSCS-01 plate. The results presented in Fig. 5(b) show a comparison of the postbuckling response for the integer-order (*α* = 1.0) CSCS-01 plate with similar results available in the literature following the classical theory of elasticity [55]. The direct comparison of the numerical results with those available in the literature for both boundary conditions shows a very good agreement. This convergence is clearly evident for load–displacement curves in the postcritical regime of the elastic response with the increasing load factor. In these results, we note that the load–displacement curves evaluated at *α* = 1.0 converge to the corresponding analytical results from classical elasticity. While this convergence is excellent for points away from the critical (turning) points, slight differences are noted for points in the neighborhood of critical points. This difference can be attributed to the sensitivity of the plate to imperfections resulting in lower values of the critical load for the imperfect structure. Unlike the current study where we consider an eccentric load *F*_{3}, Yamaki [55] studied the bifurcation of a perfect structure subject to in-plane loads. Therefore, the case considered in the present study of a critical load of an isotropic plate subject to eccentric load is lower when compared to that of the perfect plate [55]. Moreover, this difference increases with an increase in the amplitude of *F*_{3}. Therefore, the lower load factor for critical (turning) points of plates subject to eccentric loads is in accordance with the theoretical predictions. Furthermore, as expected, this difference caused by imperfections vanishes for load–displacement curves away from the critical point. In conclusion, we establish a successful convergence of the numerical method and can proceed to use it for the study of the postbuckling response of fractional-order isotropic plates.

### 4.2 Postbuckling Response of Fractional-Order Plates.

In this section, we study the effect of nonlocal interactions on the postbuckling response of fractional-order plates. Specifically, we intend to analyze the effect of the fractional-order constitutive parameters: fractional-order *α*, and nonlocal length scale *l*_{f} on (1) the critical point (turning point or bifurcation point) and (2) on both the postbuckling response and the postcritical stiffness of the structure. This analysis provides an understanding of the effects of the nonlocal behavior on the postbuckling response of thin plates. Note that either an increase in the nonlocal length scale *l*_{f} or a reduction of the fractional-order *α* correspond to higher degree of nonlocality.

The nonlinear load–displacement curves of the isotropic nonlocal plate subject to uniaxial compression *F*_{1} and an eccentric load *F*_{3} for different values of the fractional-order *α* are compared in Fig. 6. The response for the SSSS-01 plate and CSCS-01 plate are presented in Figs. 6(a) and 6(b), respectively. These figures underline a softening response of the plate, beyond the critical point (turning point), when subject to an increasing degree of nonlocality (i.e., a progressive reduction of the fraction-order *α*). This observation is valid for both types of boundary conditions being considered. This observed softening in the postbuckling nonlinear response is consistent with the behavior expected from integral-type nonlocal constitutive relations [14], particularly the fractional-order continuum theories employed [39–41].

Unlike the previous comment on the effect of the nonlocality on the postbuckling response, contrasting observations can be drawn when focusing on the effect on the critical load (bifurcation point of the nonlinear load–displacement curves). In the case of the SSSS-01 beam depicted in Fig. 6(a), a marginal increase in the critical load is observed following the increasing degree of nonlocality. At the same time, a reduction in the critical load is noted with the reducing fractional-order *α* for the CSCS-01 plate. These observations for the two case studies are in complete agreement with the study on linear buckling of fractional-order structures [49]. The fractional-order strain–displacement relations introduce fractional-order derivatives into the expressions of both the material and the geometric stiffness. This is clear from Eq. (10) for linear buckling where the fractional-order nonlocal terms $e~$ and $q~$ are present in both the numerator (material stiffness) and the denominator (geometric stiffness). An increase in the degree of nonlocality reduces the numerical values of the fractional derivatives. Thus, an increasing degree of nonlocality would result in lower material and geometric stiffness terms of the nonlocal structure. While a reduction in the material stiffness would reduce the critical load, a similar reduction in the geometric stiffness leads to higher values of the critical load. A comprehensive account of these effects, including detailed parametric studies on the effect of each stiffness term, is available in Ref. [49]. Thus, we conclude that the effect of the nonlocal interactions on the critical load is dictated by the net result of the contrasting effects caused by the reduction of individual stiffness terms. This observation, drawn in the context of the *linear* buckling analysis of nonlocal structures, is also applicable to the current discussion on *nonlinear* buckling. Based on the eigenproblem for nonlinear buckling of fractional-order plates given in Eq. (9), similar conclusions may be drawn over the effect of fractional-order derivatives on the nonlinear material and geometric stiffness terms. Therefore, the effect of the nonlocal interactions on the nonlinear buckling load depends on the net result of the fractional-order terms within the two nonlinear stiffness terms.

It is known that the slope of the load–displacement curves indicate the stiffness of the structure [2]. On the basis of the results presented in Fig. 6, we note a significant reduction in the slope of the nonlinear load–displacement curves for response beyond the critical point. This observation for *α* = 1.0 is corroborated by the existing literature on postbuckling behavior based on the classical theory of elasticity. In case of fractional values for *α*, the slope of the load–displacement curves corresponding to nonlocal elasticity is “flatter” when compared with the local analogs. This behavior confirms the lower stiffness value of nonlocal structures beyond the critical point. From the postbuckling nonlinear load–displacement curves, we also note that the influence of nonlocal elasticity increases by moving farther into the postbuckling regime. This is an important observation that points toward reduced resistance to deformation for nonlocal structures in a deep postbuckling regime.

The effect of the nonlocal length scale *l*_{f} on the postbuckling response of fractional-order isotropic plates is explored in Fig. 7. This behavior is shown in Fig. 7(a) for the SSSS-01 plate and in Fig. 7(b) for the CSCS-01 plate. Similar to the aforementioned discussion regarding the effect of the fractional-order *α*, a softening behavior in the postbuckled regime is noted for increasing the values of the nonlocal length scale, *l*_{f}. This observation is consistent for both boundary conditions chosen for this study. Furthermore, in agreement with the aforementioned results, contrasting observations are noted in terms of the effect of nonlocality on the critical load. For the SSSS-01 plate (Fig. 7(a)), the critical load increases with an increasing degree of nonlocality (higher *l*_{f}). However, the critical load of the CSCS-01 plate (Fig. 7(b)) reduces for higher values of the length scale *l*_{f}. These contrasting observations are consistent with the previous discussion on the influence of the fractional-order *α* on the critical load, as well as with previous work on linear buckling load of fractional-order plates [49].

Finally, we intend to undertake an analysis of the effect of nonlocal interactions on the sensitivity to imperfections of the plates. A detailed analysis can be proposed following the Koiter’s asymptotic approach previously carried out in Sec. 2.3 by including additional deformation energy terms corresponding to (geometric or load) imperfections in Eq. (7). However, for the sake of brevity, we present a quantitative analysis through the case study considered above. As discussed earlier, the transverse load *F*_{3} acting on the top surface of the plate amounts to a load imperfection on the isotropic plate subject to uniaxial compressive loads *F*_{1} applied the edges of the plate. Note that the role of the eccentric force *F*_{3} is to activate the stable branch of the load–displacement equilibrium curve beyond the critical point. This is because the fundamental curve *C*_{0} is no longer stable for equilibrium points beyond the critical point (see Sec. 2.3). Another approach in terms of geometric imperfection can be similarly envisioned. The distribution of *F*_{3} is motivated by the requirement to activate the transverse displacement mode corresponding to primary mode of buckling [49]. A consequence of the imperfection would be to reduce the critical load corresponding to nonlinear buckling, commonly referred to as the imperfection sensitivity [2,56]. To compare this effect for nonlocal structures, we present the load–displacement curves of the CSCS-01 fractional-order plate for different amplitudes of the transverse force $F\xaf3$ in Fig. 8. Irrespective of the amplitude of the transverse force, the postbuckling curves converge for load factor Λ away from the critical point. However, a lower nondimensionalized critical load (turning point) is noted for the increasing amplitude of *F*_{3} over the plate. This points to the imperfection sensitivity of the nonlocal plate. In this figure, we also compare the effect of nonlocal elasticity over the imperfection sensitivity for two different sets of fractional-order constitutive parameters. Any perceptible change in the critical load is not evident from the comparison between Figs. 8(a) and 8(b) for different choices of fractional-order and nonlocal length scales. Therefore, we conclude that postbuckling response of nonlocal plates exhibits imperfection sensitivity, which, however, is not influenced by the nonlocal interactions.

## 5 Conclusions

This article presented the theoretical approach and a numerical methodology for the analysis of the postbuckling elastic response of fractional-order structures. The approach is enabled by the thermodynamically consistent and positive semi-definite form of the deformation energy density that can be achieved via the fractional-order continuum theory for nonlocal structures. Further, geometrically nonlinear models of fractional-order strain–displacement relations are integrated in the formulation to enable the stability and postcritical response analyses. This theoretical approach opens the way to using energy methods to identify critical points of nonlinear equilibrium. In addition, it allows considering the influence of nonlocal interactions on both the material and the geometric stiffness terms. The possibility of using energy-based methods, afforded by the fractional formulation, allowed the application of Koiter’s asymptotic approach to ascertain the nature of the singularities exhibited by the nonlinear and nonlocal elastic response. Based on the fractional-order nonlinear stability theory, quantitative and qualitative studies were performed to understand the response of fractional-order (nonlocal) Kirchhoff plates in a postbuckling regime. The fractional-order finite element method (f-FEM) was extended to include an incremental-iterative solver and a path-following method to solve the nonlinear system of algebraic equations necessary to obtain the nonlinear equilibrium load–displacement curves. The effect of the nonlocal interactions on the nonlinear critical load is not trivial due to a simultaneous reduction of the material and geometric stiffness terms caused by the fractional-order kinematic relations. A reduction of the postbuckling stiffness due to the nonlocal interactions is noted, and it is associated with a softening response when compared with the local elastic analogs. Finally, it is established in literature that the frame-invariant fractional-order continuum theory serves as a mathematically, physically, and thermodynamically consistent alternative to the existing integer-order theories of nonlocal elasticity. Moreover, energy methods available for the fractional-order theory make it suitable to developing models that provide a better account of the scale-dependent effects on the stability and postbuckling response of complex structures. Following this, the current study opens an opportunity to perform stability of sandwich composites, layered and porous media, and biological materials like tissues and bones.

## Acknowledgment

The following work was supported by the Defense Advanced Research Project Agency (DARPA) (Grant No. D19AP00052) and the National Science Foundation (NSF) (Grant Nos. DCSD #1825837 and MOMS #1761423). The content and information presented in this manuscript do not necessarily reflect the position or the policy of the government. The material is approved for public release; distribution is unlimited.

## Conflict of Interest

There are no conflicts of interest.