    # Fractional-Order Models for the Static and Dynamic Analysis of Nonlocal Plates

This study presents the analytical formulation and the finite element solution of fractional order nonlocal plates under both Mindlin and Kirchoff formulations. By employing consistent definitions for fractional-order kinematic relations, the governing equations and the associated boundary conditions are derived based on variational principles. Remarkably, the fractional-order nonlocal model gives rise to a self-adjoint and positive-definite system that accepts a unique solution. Further, owing to the difficulty in obtaining analytical solutions to this fractional-order differ-integral problem, a 2D finite element model for the fractional-order governing equations is presented. Following a thorough validation with benchmark problems, the 2D fractional finite element model is used to study the static as well as the free dynamic response of fractional-order plates subject to various loading and boundary conditions. It is established that the fractional-order nonlocality leads to a reduction in the stiffness of the plate structure thereby increasing the displacements and reducing the natural frequency of vibration of the plates. Further, it is seen that the effect of nonlocality is stronger on the higher modes of vibration when compared to the fundamental mode. These effects of the fractional-order nonlocality are noted irrespective of the nature of the boundary conditions. More specifically, the fractional-order model of nonlocal plates is free from boundary effects that lead to paradoxical predictions such as hardening and absence of nonlocal effects in classical integral approaches to nonlocal elasticity. This consistency in the predictions is a result of the well-posed nature of the fractional-order governing equations that accept a unique solution.

## Authors

01/19/2020

### A Ritz-based Finite Element Method for a Fractional-Order Boundary Value Problem of Nonlocal Elasticity

We present the analytical formulation and the finite element solution of...
06/09/2020

### A refined dynamic finite-strain shell theory for incompressible hyperelastic materials: equations and two-dimensional shell virtual work principle

Based on previous work for the static problem, in this paper we first de...
03/03/2020

### Thermoelastic Response of Fractional-Order Nonlocal and Geometrically Nonlinear Beams

This study presents both the theoretical framework and the finite elemen...
08/26/2020

### Variable-Order Approach to Nonlocal Elasticity: Theoretical Formulation and Order Identification via Deep Learning Techniques

This study presents the application of variable-order (VO) fractional ca...
12/24/2021

### Multiscale Nonlocal Elasticity: A Distributed Order Fractional Formulation

This study presents a generalized multiscale nonlocal elasticity theory ...
08/15/2020

### Fractional-Order Structural Stability: Formulation and Application to the Critical Load of Slender Structures

This study presents the framework to perform a stability analysis of non...
04/08/2021

### Displacement-Driven Approach to Nonlocal Elasticity

This study presents a physically consistent displacement-driven reformul...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

In the recent years, following the rapid growth in many engineering fields including, but not limited to, functionally graded materials (FGM), composites, nanotechnology, and MEMS, the modeling of the static as well as the dynamic response of complex slender structures has received considerable attention. FGMs udupa2014functionally; naebe2016functionally and sandwiched designs have found many useful applications in the design of macroscale structures such as those involved in naval and automotive systems, as well as lightweight structures, such as those employed in space and aeronautic applications. Similarly, micro- and nano-structures such as thin films, carbon nanotubes, monolayer graphene sheets, and micro tubules have demonstrated far-reaching applications in atomic devices, micro/nano-electromechanical devices, sensors, and even biological implants wang2011mechanisms; pradhan2009small; narendar2011buckling. Many of these complex structures often employ a combination of slender and thin-walled structures like beams, plates, and shells, whose response was often shown to be significantly impacted by size-dependent effects, also referred to as nonlocal effects. In the case of macrostructures, these nonlocal effects were shown to be originated from material heterogeneities, interactions between layers (e.g. FGMs or composites) or unit cells (e.g. periodic media) romanoff2016using; tarasov2013review; hollkamp2019analysis; patnaik2019generalized. In other terms, nonlocal governing equations for macrostructures often result from a process of homogenization of the initial inhomogeneous system. In the case of nano- and micro-structures, these size-dependent effects have been traced back to the existence of surface and interface stresses due to either nonlocal atomic or Van der Waals interactions wang2011mechanisms. Further, geometric effects such as, for example, changes in curvature have also been shown to induce nonlocal size-dependent effects sudak2003column; pradhan2009small; sahmani2018nonlocal. It appears that the ability to accurately model these nonlocal effects as part of the structural response is paramount in many engineering applications. In order to capture scale effects, nonlocal continuum theories were developed.

From a general perspective, nonlocal continuum theories enrich the classical (local) governing equations with information of the behavior of points within a prescribed distance; the latter typically indicated as the horizon of influence or horizon of nonlocality. The key principle behind nonlocal theories relies on the idea that all the particles located inside the horizon influence one another by means of long range cohesive forces. Seminal works from Kron̈er kroner1967elasticity, Eringen eringen1972nonlocal, and several other authors nowinski1984nonlocal; lim2015higher have explored the role of nonlocality in elasticity and laid its theoretical foundation. The mathematical description of nonlocal continuum theories relies on the introduction of additional contributions in terms of gradient or integrals of the strain field in the constitutive equations. This leads to so-called “weak” gradient methods or “strong” integral methods, respectively. Gradient elasticity theories peerlings2001critical; sidhardh2018exact; sidhardh2018element account for the nonlocal behavior by introducing strain gradient dependent terms in the stress-strain constitutive law. Integral methods polizzotto2001nonlocal; bavzant2002nonlocal model nonlocal effects by defining the constitutive law in the form of a convolution integral between the strain and the spatially dependent elastic properties over the horizon of nonlocality.

Several researchers have used the above mentioned nonlocal theories to model the response of nonlocal beams and plates with particular attention to their applications in micro- and nano-devices lu2007non; reddy2010nonlocal where nonlocal effects are typically more noticeable. In particular, the effect of the nonlocality on the buckling load of the plates pradhan2009small; narendar2011buckling; sahmani2018unified as well as the vibration response of the plates barati2018nonlinear have been extensively studied. In this context, Challamel et alchallamel2016buckling; zhang2014eringen have developed a phenomenological nonlocal model from a lattice-based nonlocal model from analytic continuation of the plate lattice equations and have used it to analyze the effect of nonlocality on the buckling loads of a nonlocal plate. All these studies have shown that the introduction of nonlocality leads to a decrease in the stiffness of the structure which translates to higher static displacements yan2015exact; sahmani2018nonlocal, lower buckling loads, and lower frequencies of vibration of the nonlocal structures. We merely note that a mixture of analytical duan2007exact; yan2015exact as well numerical methods have been employed to determine the response of the nonlocal plates in the aforementioned studies. The numerical strategies included perturbation methods barati2018nonlinear; sahmani2018nonlocal, differential quadrature methods pradhan2009small; golmakani2014nonlinear, Galerkin methods phadikar2010variational; anjomshoa2013application as well as spectral collocation methods ma2016vibration.

Although these classical studies on nonlocal elasticity have been able to address several features of the response of size-dependent nonlocal structures, they encounter some key shortcomings. Gradient theories provide a satisfactory description of the material micro structure, but they introduce serious difficulties when enforcing the boundary conditions associated with the strain gradient-dependent terms peerlings2001critical; aifantis2003update

. On the other side, the integral methods are better suited to deal with boundary conditions but require the attenuation functions to have a positive Fourier transform everywhere in order to avoid instabilities

bavzant1984instability; bavzant2002nonlocal. Additionally, in both these classes of methods, the stress at any point cannot be obtained unless the strains in the neighbourhood of the particular point are known. In other terms, there exists no explicit relation for obtaining the stress at a given point from the strain at that particular point. This prevents the application of variational principles in these theories. More specifically, the basis of variational formulation is the principle of minimum total potential energy which is valid under the assumption that the stress at a point can be uniquely defined in terms of the strain at that point. Since an explicit relation between stress and strain components at a reference point cannot be found in the classical nonlocal theories, the principle of minimum potential energy cannot be applied to these theories. This aspect prevents the development of variational finite element methods to obtain numerical solutions to the classical nonlocal formulations. We note that the development of numerical methods is essential in this class of problems because the complex nonlocal governing equations associated with both gradient and integral formulations do not generally admit closed-form analytical solutions. Several researchers have developed inverse approaches in order to define a quadratic form of the total potential energy and then use it to obtain the response of the nonlocal structures via Galerkin or Ritz approximations phadikar2010variational; anjomshoa2013application. In addition to the above shortcomings, classical integral approaches lead to mathematically ill-posed governing equations which leads to erroneous predictions such as the absence of nonlocal effects and the occurrence of hardening behavior for certain combinations of boundary conditions romano2017constitutive; romano2018nonlocal. In this class of problems, the ill-posedness stems from the fact that the constitutive relation between the bending field and the curvature is a Fredholm integral of the first kind, whose solution does not generally exists and, if it exists, it is not necessarily unique romano2017constitutive; romano2018nonlocal.

In recent years, fractional calculus has emerged as a powerful mathematical tool to model a variety of nonlocal and multiscale phenomena. Fractional derivatives, which are a differ-integral class of operators, are intrinsically multiscale and provide a natural way to account for nonlocal effects. In fact, the order of the fractional derivative dictates the shape of the power-law influence function (kernel of the fractional derivative) while its interval defines the horizon of its influence, i.e., the distance beyond which information is no longer accounted for in the derivative. As a result, time-fractional operators enable memory effects (i.e. the response of a system is a function of its past history) while space-fractional operators can account for nonlocal and scale effects. Given the multiscale nature of fractional operators, fractional calculus has found wide-spread applications in nonlocal elasticity. Riesz-type fractional derivatives have been shown to emerge as the continuum limit of discrete systems (e.g. such as chains and lattices) with power-law long-range interactions tarasov2013review. Space-fractional derivatives have been used to formulate nonlocal constitutive laws drapaca2012fractional; carpinteri2014nonlocal; sumelka2014fractional; sumelka2015non; hollkamp2020application as well as to account for microscopic interaction forces cottone2009elastic; di2008long. Space-fractional derivatives have also been employed to capture attenuation including a variety of conditions such as interatomic nonlocal forces cottone2009elastic, nonlocal stress-strain constitutive relations atanackovic2009generalized, and even bandgaps in periodic media hollkamp2019analysis. Previous works conducted on the development of nonlocal continuum theories based on fractional calculus have highlighted that the differ-integral nature of the fractional operators allows them to combine the strengths of both gradient and integral based methods while at the same time addressing a few important shortcomings of the integer order formulations drapaca2012fractional; carpinteri2014nonlocal; sumelka2014fractional; patnaik2019generalized; patnaik2019FEM; sidhardh2020geometrically.

In this study, we build upon the fractional-order nonlocal continuum model proposed in patnaik2019generalized to develop a fractional-order model for nonlocal plates. The overall goal of this study is three fold. First, we derive the governing equations for the nonlocal plates in a strong form by using variational principles. We will show that the fractional-order formulation allows the application of variational principles because the stress in the fractional-order formulation can be uniquely and explicitly determined from the strain at the particular point. In fact, nonlocality will be introduced into the plate through nonlocal fractional-order kinematics and not through the classical integral constitutive relations between the stress and the strain. Second, we formulate a fully consistent and highly accurate 2D fractional-order finite element method (f-FEM) by extending the work in patnaik2019FEM on 1D fractional-order beams, to numerically investigate the response of the fractional-order nonlocal plates. Although several FE formulations for fractional-order equations have been proposed in the literature, they are based on Galerkin or Petrov-Galerkin methods that are capable of solving 1D hyperbolic and parabolic differential equations involving transport processes jin2016petrov; li2017galerkin. We develop a Ritz FEM that is capable of obtaining the numerical solution of the fractional-order plate governing equations. We highlight here that although Ritz FEMs for classical nonlocal elasticity problems have been developed in the literature, they do not extend to fractional-order nonlocal modeling, because the attenuation function capturing the nonlocal interactions in the fractional-order model involves a singularity within the kernel podlubny1998fractional. Moreover, the necessity and flexibility of the f-FEM developed here becomes clear from the complexities involved in the handling of integral boundary conditions fernandez2016bending. Finally, we use the developed 2D f-FEM to analyze the effect of nonlocality on the static and free vibration response of the fractional-order plates. We will show that independently from the boundary conditions, the fractional-order theory predicts a consistent softening behaviour for the fractional-order plates as the nonlocality degree increases. This latter aspect is unlike the paradoxical predictions of hardening or of the absence of nonlocal effects predicted by classical integral nonlocal approaches challamel2014nonconservativeness; khodabakhshi2015unified; romano2017constitutive; romano2018nonlocal; barretta2019stress for certain combinations of boundary conditions.

The remainder of the paper is structured as follows: first, we introduce the fractional-order formulation used in this study to model nonlocal elasticity. Next, we derive the governing equations of fractional-order Mindlin as well as fractional-order Kirchoff plates in strong form using variational principles. Further, we derive a strategy for obtaining the numerical solution to the plate governing equations using 2D f-FEM. Finally we validate the 2D f-FEM, establish its convergence, and then use it to analyze the effect of the fractional-order nonlocality on the static and free vibration response of plates under different types of loading conditions.

## Ii Nonlocal Elasticity via Fractional Calculus

Previous works conducted on the development of nonlocal continuum theories based on fractional calculus have highlighted its ability to combine the strengths of both gradient and integral based methods while simultaneously addressing a few important shortcomings of these integer-order formulations drapaca2012fractional; carpinteri2014nonlocal; sumelka2014fractional; patnaik2019generalized; patnaik2019FEM; sidhardh2020geometrically. Recall that the key shortcomings included the requirement of higher-order essential boundary conditions in gradient based methods, the need for a kernel that is always positive in integral methods, and the inability of both methods to leverage variational principles. To this regard, note that the kernel used in fractional derivatives is positive everywhere podlubny1998fractional. Unlike gradient elasticity methods, additional essential boundary conditions are not required when using Caputo fractional derivatives hollkamp2019analysis; patnaik2019generalized

. Further, we will show how variational principles can be used in the fractional-order formulation of nonlocality by using a specific strategy that involves defining a fractional-order (nonlocal) deformation gradient tensor. The nonlocal plate theory presented in this work builds on the fractional-order nonlocal continuum formulation presented in

patnaik2019generalized. In the following, we briefly review the fractional-order nonlocal continuum model presented in patnaik2019generalized which provides the foundation to develop the fractional-order nonlocal plate theory.

In analogy with the traditional approach to continuum mechanics, we perform the deformation analysis of a nonlocal solid by introducing two configurations, namely, the reference (undeformed) and the current (deformed) configurations. The motion of the body from the reference configuration (denoted as X) to the current configuration (denoted as x) is assumed as:

 x=Ψ(X,t) (1)

such that is a bijective mapping operation. The relative position of two point particles located at and in the reference configuration of the nonlocal medium is denoted by (see Fig. (1)). After deformation due to motion , the particles occupy the new positions and

, such that the relative position vector between them is

. It appears that and are the material and spatial differential line elements in the nonlocal medium, conceptually analogous to the classical differential line elements and in a local medium representation. Figure 1: (a) Schematic indicating the infinitesimal material and spatial line elements in the nonlocal medium under the displacement field u. (b) Horizon of nonlocality and length scales at different material points. The nonlocal model can account for a partial (i.e. asymmetric) horizon condition that occurs for points X close to a boundary or an interface. The schematic is adapted from patnaik2019generalized.

The mapping operation described in Eq. (1) captures the nonlocal behavior of the solid by leveraging fractional calculus formulation. More specifically, the differential line elements of the nonlocal medium are modeled by imposing a fractional-order transformation on the classical differential line elements as follows:

 d~x=[DαXΨ(X% ,t)]dX=[~FX(X,t)]dX (2a) d~X=[DαxΨ−1(x,t)]dx=[~Fx(%x,t)]dx (2b)

where is a space-fractional derivative whose details will be presented below. Given the differ-integral nature of the space-fractional derivative, the differential line elements and have a nonlocal character. Using the definitions for and , the fractional deformation gradient tensor with respect to the nonlocal coordinates is obtained in patnaik2019generalized as:

 d~xd~X=αF=~FXF−1~F−1x (3)

where F is the classical deformation gradient tensor given as , in local and integer-order form.

The space-fractional derivative is taken according to a Riesz-Caputo (RC) definition with order defined on the interval and given by:

 (4a) DαXjΨi(X,t)=12Γ(2−α)[Lα−1Aj CXAjDαXjΨi(X,t)−Lα−1Bj CXjDαXBjΨi(%X,t)] (4b)

where is the Gamma function, and and are the left- and right-handed Caputo derivatives of , respectively. Before proceeding, we discuss certain implications of the above definition of the fractional-order derivative. The interval of the fractional derivative defines the horizon of nonlocality (also called attenuation range in classical nonlocal elasticity) which is schematically shown in Fig. (1) for a generic point . The length scale parameters and ensure the dimensional consistency of the deformation gradient tensor, and along with the term ensure the frame invariance of the constitutive relations patnaik2019generalized. As discussed in patnaik2019generalized, the deformation gradient tensor introduced via Eq. (3) enables an efficient and accurate treatment of the frame invariance in presence of asymmetric horizons, material boundaries, and interfaces (see Fig. (1)).

In analogy with the classical strain measures, the nonlocal strain can be defined using the fractional-order differential line elements as . Using Eq. (3), the Lagrangian strain tensor in the nonlocal medium is obtained as:

 αE=12(αFTαF−I) (5)

where I is the identity tensor. Using kinematic position-displacement relations, the expressions of the strains can be obtained in terms of the displacement gradients. The fractional displacement gradient is obtained using the definition of the fractional deformation gradient tensor from Eq. (1) and the displacement field as:

 ∇αUX=~FX−I (6)

The fractional gradient denoted by is given as . Using the nonlocal strain defined in Eq. (5) and the fractional deformation gradient tensor given in Eq. (3) together with Eq. (6), the relationship between the infinitesimal strain tensor and displacement gradient tensor is obtained as:

 ϵ=12(∇αUX+∇αUTX) (7)

For small displacements, the above infinitesimal strain tensor has the same definition in both Eulerian and Lagrangian descriptions patnaik2019generalized. Further, the stress tensor in the fractional-order model of the nonlocal medium has a form analogous to the local case as:

 σij=Cijklϵkl (8)

where denotes the constitutive matrix of the solid. We emphasize that the stress defined through the above equation is nonlocal in nature. This follows from the fractional-order definition of the deformation gradient tensor which is then reflected in the nonlocal strain as evident from Eq. (5). As expected, classical continuum mechanics relations are recovered when the order of the fractional derivative is set as .

Note that in the above presented fractional-order formulation, nonlocality has been modeled using fractional-order kinematic relations. More specifically, differential line elements in the undeformed and deformed nonlocal configurations were modeled using fractional-order deformation gradients which, in turn, were used to obtain the strain in the nonlocal medium. This definition of the strain has critical implications on the nonlocal formulation. Note that, given the strain at a point, the stress at the same point can be uniquely and explicitly determined by using Eq. (8). Recall that the basis of variational formulation is the principle of minimum total potential energy which is valid under the assumption that the stress at a point can be uniquely defined in terms of the strain at that point. It is immediate to deduce that the fractional-order model of nonlocality allows the application of variational principles. It is exactly this opportunity offered by the fractional operators that forms the basis for the development of a fractional-order model of nonlocal plates via the application of the extended Hamilton’s principle. Additionally, this fractional-order formulation of nonlocality also ensures a quadratic form of the potential energy of the system as we will show in §III. In the following, we will use this fractional-order formulation to model the response of both nonlocal Mindlin and Kirchoff plates.

## Iii Fractional-Order Modeling of Nonlocal Mindlin Plates

We use the fractional-order continuum formulation presented above to develop a fractional-order analogue of the classical Mindlin plate formulation. A schematic of the undeformed rectangular plate along with the chosen Cartesian reference frame is illustrated in Fig. (2). The top surface of the plate is identified as , while the bottom surface is identified as . The domain corresponding to the mid-plane of the plate (i.e., ) is denoted as , such that where and are the length and width of the plate, respectively. The domain of the plate is identified by the tensor product . The edges forming the boundary of the mid-plane of the plate are denoted as as shown in the Fig. (2). Figure 2: Schematic of the rectangular plate showing key geometric parameters. The midplane of the plate is indicated by the z=0 plane of the domain (Ω). The boundaries of the domain are also indicated in the schematic.

For the chosen coordinate system, the in-plane and transverse components of the displacement field, denoted by , and at any spatial location , are related to the mid-plane displacements in the following manner:

 u(x,y,z,t)=u0(x,y,t)−zθx(x,y,t) (9a) v(x,y,z,t)=v0(x,y,t)−zθy(x,y,t) (9b) w(x,y,z,t)=w0(x,y,t) (9c)

where , , and are the mid-plane displacements of the plate along the , , and directions. and are the rotations of the transverse normal about the and axes, respectively. In the interest of a more compact notation, the functional dependence of the displacement fields on the spatial and the temporal variables will be implied unless explicitly expressed to be constant. Based on the above displacement fields, the strain components in a fractional-order Mindlin plate are evaluated using Eq. (7) as:

 ϵxx=Dαxu0−zDαxθx (10a) ϵyy=Dαyv0−zDαyθy (10b) γxy=Dαyu0+Dαxv0−z(Dαyθx+Dαxθy) (10c) γxz=Dαxw0−θx (10d) γyz=Dαyw0−θy (10e) ϵzz=0 (10f)

Note that similarly to the classical Mindlin theory, the in-plane strain components vary in a linear fashion through the plate thickness, while the transverse shear strains are constant through the thickness. The corresponding stresses are determined using the linear stress-strain relationships given in Eq. (8).

By using the above defined strain and stress fields, we derive the strong-form of the governing equations for the fractional-order plate using the extended Hamilton’s principle:

 ∫t2t1(δU−δV−δT)dt=0 (11)

The nonlocal virtual strain energy , the virtual work done by externally applied forces , and the virtual kinetic energy are obtained as:

 (12a) δV=∫Ω[Fxδu0+Fyδv0+Fzδw0+Mθxδθx+Mθyδθy]dΩ (12b) (12c)

Note that for a rectangular plate. are the external loads applied in the , and directions, respectively, and

are the external moments applied about the

and axes, respectively. By substituting the expression for the stress given in Eq. (8) in Eq. (12a), it is immediate that the potential energy of the fractional-order nonlocal plate is quadratic in nature. In fact, the use of the fractional-order formulation leads to a self-adjoint and positive definite system. The proof for the same can be found in patnaik2019FEM where the fractional-order formulation introduced earlier in §II was used to model nonlocal beams. The same proof directly extends to the fractional-order plate formulations presented in the following and hence we do not provide it again here.

We first simplify the virtual strain energy in order to relieve the variations of any differentiation. By substituting the expressions for the strains in the expression for we obtain:

 δU=∫Ω[NxxδDαxu0+MxxδDαxθx+NyyδDαyv0+MyyδDαyθy+Nxyδ(Dαyu0+Dαxv0)+Mxyδ(Dαyθx+Dαxθy)+Qxzδ(Dαxw0−θx)+Qyzδ(Dαyw0−θy)]dΩ (13)

where the in-plane stress resultants , the transverse stress resultants and the moment resultants are defined as:

 (14a) {Mxx,Myy,Mxy}=∫h2−h2{−zσxx,−zσyy,−zσxy}dz (14b)

where is the shear correction factor. The expression for in Eq. (13) is further simplified using integration by parts and the definitions of the fractional derivatives. Here below, we provide the result of this simplification for only one term within the integral in Eq. (13):

 ∫ΩNxx[δDαxu0]dxdy=−∫Ω[DαxNxx]δu0dΩ+∫Γx[I1−αxNxx]δu0dy (15)

The detailed steps leading to the above simplification can be found in patnaik2019FEM where a similar variational approach has been used in the context of 1D beams. In the above Eq. (15), is the Riesz Riemann-Liouville derivative of order which is defined as:

 Dαxψ=12Γ(2−α)[lα−1Bx(RLx−lBxDαxψ)−lα−1Ax(RLxDαx+lAxψ)] (16)

where is an arbitrary function and and are the left- and right-handed Riemann Liouville derivatives of to the order , respectively. The Riesz fractional integral is defined as:

 I1−αxψ=12Γ(2−α)[lα−1Bxx−lBxI1−αxψ−lα−1AxxI1−αx+lAxψ] (17)

where and are the left and right Riesz integrals to the order , respectively. Note that the fractional derivative and the fractional integral are defined over the interval unlike the fractional derivative which is defined over the interval . This change in the terminals of the interval of the Riesz Riemann-Liouville fractional integral and derivative follows from the integration by parts technique used to simplify the variational integrals (see patnaik2019FEM). The remaining terms within the integral in Eq. (13) can be simplified in a similar fashion.

Note that the expression for the virtual work given in Eq. (12b) can be directly used within Eq. (11). Further, the expression of the virtual kinetic energy matches exactly its counterparts in the classical integer-order plate theory (see reddy2006theory) and can be expressed as:

 δT=−∫Ω[I0¨u0δu0+I0¨v0δv0+I0¨w0δw0+I2¨θxδθx+I2¨θyδθy]dΩ (18)

where and .

The expressions for the virtual quantities are substituted within the Hamilton’s principle statement in Eq. (11). The governing equations and the boundary conditions of the fractional-order Mindlin plate are then obtained by using the fundamental lemma of variational calculus as:

 DαxNxx+DαyNxy+Fx=I0∂2u0∂t2 (19a) DαxNxy+DαyNyy+Fy=I0∂2v0∂t2 (19b) DαxQxz+DαyQyz+Fz=I0∂2w0∂t2 (19c) DαxMxx+DαyMxy−Qxz+Mθx=I2∂2θx∂t2 (19d) DαxMxy+DαyMyy−Qyz+Mθy=I2∂2θy∂t2 (19e)

The corresponding essential and natural boundary conditions are obtained as:

 δu0=0, δv0=0, δw0=0, δθx=0, δθy=0  ∀ Γx∪Γy (20a) I1−αxNxx=0, I1−αxNxy=0, I1−αxQxz=0, I1−αxMxx=0, I1−αxMxy=0  ∀ Γx (20b) I1−αyNxy=0, I1−αyNyy=0, I1−αyQyz=0, I1−αyMxy=0, I1−αyMyy=0  ∀ Γy (20c)

Note that the natural boundary conditions are nonlocal in nature. This is similar to what is seen in classical integral approaches to the modeling of nonlocal plates lu2007non. We anticipate that the nonlocal nature of the natural boundary conditions does not concern us immediately as we will solve the above system of equations using a FE technique. Recall that natural boundary conditions are implicitly satisfied when obtaining the solutions using FE techniques and are accurate up to the order of the specific FEM. Additionally, the following initial conditions are required to obtain the transient response:

 δu0=0, δv0=0, δw0=0, δθx=0, δθy=0  ∀ Ω at t=0 (21a) δ˙u0=0, δ˙v0=0, δ˙w0=0, δ˙θx=0, δ˙θy=0  ∀ Ω % at t=0 (21b)

Note that the governing equations for the in-plane and transverse displacements are uncoupled, similar to what is seen in the classical Mindlin plate formulation. As expected, the classical plate governing equations and boundary conditions are recovered for . Note that the solution of the above equations yield the mid-plane displacements. The entire displacement field of the plate can then be obtained using Eq. (9).

The plate governing equations given in Eq. (19) can be expressed in terms of the displacement field variables by using the constitutive stress-strain relations of the plate. For the sake of generality, here below, we provide the expressions for an orthotropic plate:

 A11Dαx[Dαxu0]+A12Dαx[Dαyv0]+A66Dαy[Dαyu0+Dαxv0]+Fx=I0∂2u0∂t2 (22a) A66Dαx[Dαyu0+Dαxv0]+A12Dαy[Dαxu0]+A22Dαy[Dαyv0]+Fy=I0∂2v0∂t2 (22b) (22c) (22d) D66Dαx[Dαyθx+Dαxθy]+D12Dαy[Dαxθx]+D22Dαy[Dαyθy]+KsA44[Dαyw0−θy]+I2∂2θy∂t2=Mθy (22e)

where the different material constants used in the above equations are given as:

 A11=E1h1−ν12ν21, A12=ν12E2h1−ν12ν21, A22=E2h1−ν12ν21, A44=G23h, A55=G13h, A66=G12h (23)

In the above equation, and are the moduli of elasticity along the and axes, respectively. and are the Poisson’s ratios. Recall that, for an orthotropic material, we have . We highlight here that the different subscripts used with the above material constants are consistent with the notations used in literature. Additionally, we also have:

 D11=E1h312(1−ν12ν21), D12=ν12E2h312(1−ν12ν21), D22=E2h312(1−ν12ν21), D66=G12h312 (24)

The expressions for the boundary conditions in terms of the displacement variables are given as:

 (25)
 (26)

The governing equations for an isotropic fractional-order plate can be obtained by setting , , and in the above equations.

## Iv Fractional-Order Modeling of Nonlocal Kirchoff Plates

In a similar way to the approach used above, we apply variational principles to develop a fractional-order analogue of the classical Kirchoff plate theory. Recall that the Kirchoff theory is applicable only to thin plates, that is when the plate in-plane characteristic dimension to thickness ratio is on the order of 50 or greater reddy2006theory. Under such conditions, the transverse shear strains and can be neglected and the rotations and , introduced previously in §III, are approximated as:

 {θx, θy}≈{∂w0∂x, ∂w0∂y} (27)

The strain-displacement relations for the fractional-order Kirchoff plate are derived by substituting the above assumptions in the strain-displacement relations of the fractional-order Mindlin plate given in Eq. (9). Using the above formalism, the fractional-order strains for the Kirchoff plate are obtained as:

 ϵxx=Dαxu0−zDαx(∂w∂x) (28a) ϵyy=Dαyv0−zDαy(∂w∂y) (28b) (28c)

The corresponding stresses are determined using the stress-strain relationships given in Eq. (8). The governing differential equations and the corresponding boundary and initial conditions of the fractional-order Kirchoff plate are derived by using variational principles as illustrated in §III for fractional-order Mindlin plates, hence we do not provide the detailed derivation here. It appears that the differential equations governing the linear in-plane and the transverse response are uncoupled similar to the Mindlin plates. Given the nature of the kinematic relations in Eqs. (9,27), it is immediate that the differential equations and the boundary conditions governing the in-plane response of both the fractional-order Mindlin and the fractional-order Kirchoff plates are identical to each other, similar to what is established in classical plate theories. Hence, we do not provide them explicitly here. However, the governing equation corresponding to the transverse response of the fractional-order Kirchoff plate is different from that of the fractional-order Mindlin plate and it is obtained as:

 (29)

The corresponding boundary conditions and initial conditions are obtained as:

 (30a) (30b) ∀ Ω at t=0, δw0=0, δ˙w0=0 (30c)

In the above equations, denotes the first-integer order derivative with respect to the spatial variables , and denotes the second integer-order derivative with respect to time. The governing equations and the corresponding boundary conditions given in Eqs. (29,30) are expressed in terms of the displacement variables for an orthotropic plate as:

 D11D1+αx[D1+αxw0]+2D12D1+αx[D1+αyw0]+D22D1+αy[D1+αyw0]+D66[D1yDαx(DαxD1yw0+DαyD1xw0)+D1xDαy(DαxD1yw0+DαyD1xw0)]=Fz−I0¨w0+I2∇2¨w0 (31a) ∀ Γx⎧⎪ ⎪⎨⎪ ⎪⎩[l]δw0=0orD11Dαx[D1+αxw0]+2D66[Dαy(DαxD1yw0+DαyD1xw0)]+D12[DαxD1+αyw0]−I2D1x¨w0=0δD1xw0=0orD11D1+αxw0+D12D1+αyw0=0 (31b) ∀ Γy⎧⎪ ⎪⎨⎪ ⎪⎩[l]δw0=0orD22Dαy[D1+αyw0]+2D66[Dαx(DαxD1yw0+DαyD1xw0)]+D12Dαy[D1+αxw0]−I2D1x¨w0=0δD1yw0=0orD12D1+αxw0+D22D1+αyw0=0 (31c)

We highlight here that we have avoided using brackets within the various operators in the above equations for the sake of brevity. We emphasize that the various operators in the above equation are applied sequentially, for example, .

## V 2D Fractional Finite Element Method (f-FEM)

The fractional-order nonlocal governing equations for the nonlocal Mindlin and Kirchoff plates are numerically solved via a nonlocal finite element method. The FE formulation developed for the fractional-order governing equations builds upon the FE methods developed in the literature for integral models of nonlocal elasticity pisano2009nonlocal; phadikar2010variational; norouzzadeh2017finite. However, several modifications are necessary owing to both the choice and the behavior of the attenuation functions used in the definition of the fractional-order derivatives, as well as to the nonlocal continuum model adopted in this study. In the following, we present the f-FEM for the fractional-order Mindlin theory and then we outline the modifications that would be required in the f-FEM for the fractional-order Kirchoff theory.

### v.1 2D f-FEM Formulation

We formulate the f-FEM starting from a discretized form of the Hamiltonian functional given in Eq. (11). For this purpose, the plate domain is uniformly discretized into disjoint four-noded quadrilateral (Q4) elements , such that , and . is the total number of discretized elements. The Cartesian coordinates of each point

are interpolated by using the

Lagrangian shape functions for Q4 elements in the following manner:

 (x,y)=({L}i{X}ei,{L}i{Y}ei) (32)

where is a row vector consisting of the shape functions, while and are column vectors consisting of the and coordinates of the nodes of the element

. The vector containing the nodal degrees of freedom of the element

is given as:

 {Uei}={{u0 v0 w0 θx θy}i,1 {u0 v0 w0 θx θy}i,2 {u0 v0 w0 θx θy}i,3 {u0 v0 w0 θx θy}i,4}T (33)

where the subscript denotes the element number and the local node number . The unknown displacement field variables at any point are evaluated by interpolating the corresponding nodal degrees of freedom of . For example, at a point can be obtained as:

 u0(x)={L(1)i 0 0 0 0 L(2)i 0 0 0 0 L(3)i 0 0 0 0 L(4)i 0 0 0 0}{Uei}≡{L(u)i(x)}{Uei} (34)

The superscript in the row vector indicates the specific displacement variable being interpolated which is in Eq. (34) and the subscript denotes the element number. The other displacement variables can be obtained using similar interpolations.

We use the discretization scheme discussed above to approximate the Hamiltonian of the system. We start by deriving the discretized form of the nonlocal virtual strain energy. The stress and moment resultants in Eq. (13) can be expressed as:

 {Nxx,Nyy,Nxy,Mxx,Myy,Mxy}T=[SB]{Dαxu0,Dαxv0,Dαyu0+Dαxv0,Dαxθx,Dαyθy,Dαyθx+Dαxθy}T (35a) {Qyz,Qxz}T=[SS]{Dαyw0−θy,Dαxw0−θx}T (35b)

where and are the constitutive matrices of the plate. The different elements of these matrices have already been presented for orthotropic materials in Eqs. (23,24). It is immediate that the approximation of the virtual strain energy requires the approximation of the different fractional-order derivatives in Eq. (35).

Consider the fractional-order derivative at the point . Using the definition of the RC derivative given in Eq. (4), is expressed as:

 Dαx[u0(x)]=12(1−α)⎡⎣lα−1Ax∫xx−lAxD1x′[u0(x′)](x−x′)α dx′+lα−1Bx∫x+lBxxD1x′[u0(x′)](x′−x)α dx′⎤⎦ (36)

where

is a dummy variable along the

axis used within the definition of the fractional-order derivative. The above expression can be recast as:

 (37a) where K(x,x′,lAx,lBx,α)={12(1−α) lAxα−1 |x−x′|−α  x′∈(x−lAx,x)12(1−α) lBxα−1 |x−x′|−α  x′∈(x,x+lBx) (37b)

is the kernel of the fractional derivative. Note that kernel is a function of the relative distance between the points and , and can be interpreted similarly to the attenuation functions used in integral models of nonlocal elasticity. Clearly, the attenuation decays as a power-law in the distance with an exponent equal to the order of the fractional derivative. Note that contains the integer-order derivative . is evaluated at in terms of the nodal displacement variables corresponding to the element , such that . Using Eq. (34), the integer-order derivative can be expressed as:

 D1x′[u0(x′)]=[Bu,x(x′)]{Uep} (38)

where the is given as:

 (39)

The subscript in indicates that the displacement variable under consideration is and the direction of the differentiation is . Using the above expression for the integer-order derivative , the fractional-order derivative in Eq. (37) is obtained as:

 Dαx[u0(x)]=∫x+lBxx−lAxK(x,x