{'fraction_non_alphanumeric': 0.03751121259072005, 'fraction_numerical': 0.012456168963548887, 'mean_word_length': 4.7459294834250905, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Geometrically Nonlinear Response of a Fractional-Order Nonlocal Model of Elasticity

This study presents the analytical and finite element formulation of a geometrically nonlinear and fractional-order nonlocal model of an Euler-Bernoulli beam. The finite nonlocal strains in the Euler-Bernoulli beam are obtained from a frame-invariant and dimensionally consistent fractional-order (nonlocal) continuum formulation. The finite fractional strain theory provides a positive definite formulation that results in a unique solution which is consistent across loading and boundary conditions. The governing equations and the corresponding boundary conditions of the geometrically nonlinear and nonlocal Euler-Bernoulli beam are obtained using variational principles. Further, a nonlinear finite element model for the fractional-order system is developed in order to achieve the numerical solution of the integro-differential nonlinear governing equations. Following a thorough validation with benchmark problems, the fractional finite element model (f-FEM) is used to study the geometrically nonlinear response of a nonlocal beam subject to various loading and boundary conditions. Although presented in the context of a 1D beam, this nonlinear f-FEM formulation can be extended to higher dimensional fractional-order boundary value problems.


I. INTRODUCTION
Several theoretical and experimental studies have shown that nonlinear and scale-dependent effects are prominent in the response of several structures ranging from sandwich composites, to layered and porous media, to random and fractal media, to media with damage and cracks, to biological materials like tissues and bones [1][2][3][4][5][6]. In all these classes of problems, the coexistence of different material spatial scales renders the response fully nonlocal [7][8][9]. In many applications involving slender structures, large and rapidly varying loads typically induce nonlinear response. These nonlinear size-dependent effects are particularly prominent in microstructures and nanostructures that have far-reaching applications in medical devices, atomic devices, micro/nano-electromechanical devices and sensors [10][11][12]. These devices are primarily made from a combination of slender structures such as beams, plates and shells. Clearly, accurate modeling of the nonlocal and nonlinear effects in the response of these structures is paramount in many engineering as well as biomedical applications. However, it appears that there is a limited amount of studies focusing on the modeling of geometrically nonlinear and nonlocal slender structures. In the following, we briefly review the main characteristics of these studies as well as discuss key limitations.
Nonlocal continuum theories use different strategies to enrich the classical governing equations, that describe the response of the medium at a point, with information on the behavior of distant points contained in a prescribed area of influence. Depending on the field in which this nonlocal description of the medium is applied, the same concept can be referred to with different terminologies such as action at a distance, long range interactions, or scale effects. The key concept at the basis of nonlocal continuum theories is the horizon of influence or horizon of nonlocality which describes areas of the medium where distant points can still influence one another via long range cohesive forces [7][8][9]. Several theories based on integral methods, gradient methods, and very recently, the peridynamic approach have been developed to capture these long range interactions and analyze their effect on the response of either media or structures. Gradient elasticity theories [13][14][15][16][17] account for the nonlocal behavior by introducing strain gradient dependent terms in the stress-strain constitutive law. Integral methods [18][19][20] define constitutive laws in the form of convolution integrals between the strain and the spatially dependent elastic properties over the horizon of nonlocality.
As emphasized earlier, several applications involving nonlocal slender structures also experience geometrical nonlinearities. Nonlocal effects on the nonlinear transverse response of shear deformable beams were modeled in [21,22] using von Kármán nonlinearity and Eringen's nonlocal elasticity. Reddy et al. [23] modeled a nonlocal plate starting from the differential constitutive relations of Eringen and the von Kármán nonlinear strains. The differential model of Eringen's nonlocal elasticity was also used in [24,25] to analyse nonlocal and nonlinear beams. Similar models for beams and plates have also been developed in [26,27]. Several studies were also carried out to simulate the multiphysics response of beams such as, for example, the nonlinear static and dynamic response of a capacitive electrostatic nanoactuator [28] or the effect of nonlocal elasticity over the electroelastic response of flexoelectric nanobeams [20].
Although the above nonlinear studies based on classical approaches have been able to capture key aspects of nonlocal elasticity, they are subject to important drawbacks. As shown in [29][30][31], the integral nonlocal methods do not allow defining a self-adjoint quadratic potential energy which results in paradoxical predictions such as hardening and absence of nonlocal effects for certain combinations of boundary conditions [30][31][32]. Also, given the complex nature of the nonlocal and nonlinear differential equations, research on the corresponding numerical methods is limited so far. It appears that only quadrature based methods [28,33,34] have been developed for the solution of nonlinear and nonlocal boundary value problems (BVPs).
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 influence function (kernel of the fractional derivative) while its domain of integration defines the horizon of influence, that is the distance beyond which information is no longer accounted for in the derivative. It follows that time fractional operators enable memory effects (i.e. the response of a system is a function of its past history) while space fractional operators account for nonlocal and scale effects. The unique set of properties of fractional operators have determined a surge of interest in the exploration of possible applications. Among the areas that have rapidly developed, we find the modeling of viscoelastic materials [35,36], transport processes in complex media [37][38][39][40], and model-order reduction of lumped parameter systems [41,42]. Given the multiscale nature of fractional operators, fractional calculus has also found several applications in nonlocal elasticity. Space-fractional derivatives have been used to formulate nonlocal constitutive laws [43,44] as well as to account for microscopic interaction forces [45][46][47]. Very recently fractional-order theories have been extended to model and analyze the static response, buckling characteristics, as well as the dynamic response of nonlocal beams [10,12,48,49]. However, all these studies focused on the linear geometric analysis of nonlocal structures.
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 these integer-order formulations which have been highlighted earlier. In this study, we address these shortcomings by developing a consistent nonlinear and fractional-order nonlocal beam model. Further, we develop a nonlinear fractional-order finite element method (f-FEM) to obtain the numerical solution of nonlinear and nonlocal BVPs. The overall goal of this study is two-fold. First, we derive the governing equations for the geometrically nonlinear and nonlocal beam in strong form using variational principles. We build upon the fractional-order nonlocal continuum model proposed in [49] to develop a fractional-order constitutive relation for a geometrically nonlinear Euler-Bernoulli beam. The model proposed in [49] was shown to be self-adjoint and positive definite and thus overcoming one of the important limitations mentioned above. Further, note that the fractional-order nonlocal continuum model used in [49] is frame-invariant, dimensionally consistent, capable of dealing with asymmetric horizons, and it admits integer-order boundary conditions. Second, we formulate a fully consistent and highly accurate nonlinear f-FEM to numerically investigate the response of the geometrically nonlinear fractional-order nonlocal beams. 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 only linear hyperbolic and parabolic differential equations [50][51][52][53]. We develop a Ritz FEM that is capable of obtaining the numerical solution of the nonlinear fractional-order elliptic boundary value problem that describes the static response of the fractional-order nonlocal and nonlinear beam. Further, by using the f-FEM we show that, independently from the boundary conditions, the fractional-order theory predicts a softening behavior for the geometrically nonlinear fractional-order beam as the degree of nonlocality increases. With these results, we explain the paradoxical predictions of hardening and the absence of nonlocal effects for certain combinations of boundary conditions, as predicted by classical integral approaches to nonlocal elasticity [30][31][32].
The remainder of the paper is structured as follows: we begin with the development of the nonlinear fractional-order beam theory based on the fractional-order von Kármán strain-displacement relations. Next, we derive the governing equations of the nonlinear fractional-order beam in strong form using variational principles. Then, we derive a strategy for obtaining the numerical solution to the nonlinear beam governing equation using fractional-order FEM. Finally, we validate the nonlinear f-FEM, establish its convergence, and then use it to analyze the effect of nonlinearity and the fractional-order nonlocality on the static response of the beam under different loading conditions.

II. NONLINEAR AND NONLOCAL FRACTIONAL CONSTITUTIVE MODEL
In this section, the geometrically nonlinear response of a nonlocal Euler-Bernoulli beam is modeled using the fractional-order nonlocal continuum formulation from [40]. Starting from this continuum theory, the nonlinear constitutive relations for the nonlocal Euler-Bernoulli beam are derived. Finally, the governing differential equations and the corresponding boundary conditions for the fractional-order beam theory are also developed using variational principles.

A. Fundamentals of the fractional-order nonlocal continuum formulation
In analogy with the traditional approach to continuum modeling, we perform the deformation analysis of a nonlocal solid by introducing two stationary 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: such that Ψ(X) is a bijective mapping operation. The relative position vector of two point particles located at P 1 and P 2 in the reference configuration of the nonlocal medium is denoted by dX (see Fig. (1)). After deformation due to the motion Ψ(X), the particles move to new positions p 1 and p 2 , such that the relative position vector between them is dx. It appears that dX and dx are the material and spatial differential line elements in the nonlocal medium, conceptually analogous to the classical differential line elements dX and dx. that occurs at points X close to a boundary or interface. The schematic is adapted from [40].
The mapping operation described in Eq. (1) is used to model the differential line elements in the nonlocal solid using fractional calculus in order to introduce nonlocality into the system [40]. More specifically, the differential line elements of the nonlocal medium are modeled by imposing a fractionalorder transformation on the classical differential line elements as follows: where D α is a space-fractional derivative whose details will be presented below. Given the differintegral nature of the space-fractional derivative, the differential line elements dX and dx have a nonlocal character. Using the definitions for dX and dx, the fractional deformation gradient tensor α F with respect to the nonlocal coordinates is obtained in [40] as: where F is the classical deformation gradient tensor given as F = dx/dX given in local and integer order form.
The space-fractional derivative D α X Ψ(X, t) is taken here according to a Riesz-Caputo (RC) definition with order α ∈ (0, 1) defined on the interval X ∈ (X A , X B ) ⊆ R 3 , and given by: where Γ(·) is the Gamma function, and C X A D α X Ψ and C X D α X B Ψ are the left-and right-handed Caputo derivatives of Ψ, respectively. Before proceeding, we highlight certain implications of the above definition of the fractional-order derivative. The interval of the fractional derivative (X A , X B ) defines the horizon of nonlocality (also called attenuation range in classical nonlocal elasticity) which is schematically shown in Fig. (1) for a generic point X ∈ R 2 . The length scale parameters L α−1 Aj and L α−1 Bj ensure the dimensional consistency of the deformation gradient tensor, and along with the term 1 2 Γ(2−α) ensure the frame invariance of the constitutive relations [40]. As discussed in [40], 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 ( Fig. (1)).
The differential line elements are used to define strain measures in the nonlocal medium. Here, we provide only the definition of strain measures in Lagrangian coordinates. The corresponding definitions in Eulerian coordinates can be found in [40]. In analogy with the classical strain measures, the fractional strain can be defined using the differential line elements as dxdx − dXdX. Using Eq. (3), the Lagrangian strain tensor in the nonlocal medium is obtained as: where I is the second-order 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 U(X) = x(X) − X as: The fractional gradient denoted by ∇ α U X is given as ∇ α U Xij = D α Xj U i . Using the nonlocal strain defined in Eq. (5) and the fractional deformation gradient tensor α F given in Eq. (3) together with Eq. (6), the relationship between the strain and displacement gradient tensors is obtained as: From the strain-displacement relations, assuming moderate rotations but small strains, the fractionalorder von-Kármán strain-displacement relations are defined as: where u 1 and u 3 denote the displacement fields in the x 1 and x 3 directions. In this study, we have analyzed the effect of the above nonlinear strain-displacement relations on the response of an isotropic solid. The stress in this isotropic solid is given as: where C ijkl is the constitutive matrix for the isotropic solid and˜ kl is the fractional-order strain.

B. Nonlinear constitutive relations
We use the fractional-order nonlinear strain-displacement formulation presented above to derive the constitutive model for geometrically nonlinear nonlocal response of the beam. A schematic of a beam subject to transverse loads is shown in Fig. (2). The Cartesian coordinates are chosen such that x 3 = ±h/2 coincide with the top and bottom surfaces of the beam, and x 1 = 0 and x 1 = L correspond to the two ends of the beam in the longitudinal direction. The origin of the reference frame is chosen at the left-end of the beam and x 3 = 0 is coincident with the mid-plane. In the above defined coordinate system, the axial and transverse components of the displacement field u(x 1 , x 3 ), denoted by u 1 (x 1 , x 3 ) and u 3 (x 1 , x 3 ) at any spatial location x(x 1 , x 3 ), are related to the mid-plane displacements according to the Euler-Bernoulli assumptions: where u 0 (x 1 ) and w 0 (x 1 ) are the axial and transverse displacements of a point on the mid-plane, respectively. For the above displacement field, the fractional-order nonlinear axial strain is obtained from the fractional-order von-Kármán strain-displacement relations defined in Eq. (8) as: where the fractional-order axial strain˜ 0 (x 1 ) and bending strainκ(x 1 ) are given as: In the above equation, D α x1 (·) denotes the fractional-order RC derivative given in Eq. (4) such that α ∈ (0, 1). x A (x A1 , 0) and x B (x B1 , 0) are the terminals of the left-and right-handed Caputo derivatives within the RC derivative. Consequently, the domain (x A , x B ) is the horizon of nonlocal interactions for x(x 1 , 0), along the mid-plane. Further, l A = x 1 − x A1 and l B = x B1 − x 1 are the length scales to the left and right hand sides of the point x(x 1 , 0) along the directionx 1 .
The axial stressσ 11 resulting from the fractional-order strain is evaluated using Eq. (9) as: where E denotes the Young's Modulus of the isotropic elastic beam. The deformation energy U accumulated in the nonlocal elastic beam due to the above stress-strain distribution is obtained as: where Ω denotes the volume occupied by the beam. The total potential energy functional, assuming that the beam is subjected to distributed forces F a (x 1 ) and F t (x 1 ) acting on the mid-plane along the axial and transverse directions, is given by: To obtain the total potential energy of the system, we have assumed no body forces to be applied.
Using the above constitutive model for the fractional-order nonlocal and nonlinear elastic beam, we can derive the governing differential equations in their strong form and the associated boundary conditions by imposing optimality conditions on the total potential energy functional Π[u(x)].

C. Governing equations
Using a variational approach, the governing differential equations and the associated boundary conditions for the fractional-order beam are: and subject to the boundary conditions: In the above, N (x 1 ) and M (x 1 ) are the axial and bending stress resultants, and are given by: In the fractional-order governing equations provided in Eq.
is a Riesz-type Riemann-Liouville (R-RL) derivative of order α ∈ (0, 1). The R-RL fractional derivative of an arbitrary function f (x 1 ) is defined analogous to the RC derivative given in Eq. (4) as: are the left-and right-handed Riemann Liouville (RL) derivatives of f (x 1 ) to the order α, respectively. The detailed steps to obtain the strong form of the governing equations by minimization of the total potential energy are reported in the Appendix.
Note that the governing equations for the axial and transverse displacements are coupled, similar to what is seen in the classical nonlinear Euler-Bernoulli beam formulation. This is unlike the linear elastic fractional-order response of a nonlocal beam, where the governing equations for the axial and transverse displacements are uncoupled [49]. Further, as expected, the classical nonlinear Euler-Bernoulli beam governing equations and boundary conditions are recovered for α = 1. As mentioned previously, the above form of the governing equations for the fractional-order beam do not admit closed-form analytical solutions for the most general loading and boundary conditions. Therefore, in the following, we develop a fractionalorder finite element method (f-FEM) to obtain the numerical solution of the nonlinear fractional-order governing equations. We will then use the f-FEM to investigate the characteristic response of a nonlinear fractional-order beam.

III. FRACTIONAL FINITE ELEMENT METHOD (F-FEM)
The fractional-order nonlocal and nonlinear boundary value problem described by Eqs. (16a-16e) are numerically solved using finite element method. The fractional-order boundary value problem builds upon the FE methods developed in the literature for integral models of nonlocal elasticity [54,55]. However, several modifications are necessary due to both the specific form of the attenuation functions used in the definition of the fractional-order derivatives, and to the nonlocal continuum model adopted in this study. Note that, although a f-FEM formulation was developed in [49], the study dealt exclusively with linear fractional-order boundary value problems, while here we present the approach in the context of a nonlinear BVP. We also highlight that, although the f-FEM is presented here for nonlinear fractionalorder beams, its formulation can be easily extended to model the nonlocal and nonlinear behavior of plates and shells.

A. f-FEM formulation
We formulate the f-FEM starting from a discretized form of the total potential energy functional Π[u(x)] given in Eq. (15). For this purpose, the beam domain Ω = [0, L] along the x 1 direction, is uniformly discretized into disjoint two-noded elements Ω e i = (x i 1 , x i+1 1 ), such that ∪ Ne i=1 Ω e i = Ω, and Ω e j ∩Ω e k = ∅ ∀ j = k. N e is the total number of discretized elements and the length of each element Ω e ∈ Ω is l e . The unknown displacement field variables u 0 (x 1 ) and w 0 (x 1 ) at any point x 1 ∈ Ω e i are evaluated by interpolating the corresponding nodal degrees of freedom for Ω e i . From the definitions given in Eq. (11), it is clear that the interpolation of the axial and transverse displacement fields following Euler-Bernoulli beam theory would require at least Lagrangian (L i , i = 1, 2) and Hermitian (H j , j = 1, 2, 3, 4) interpolation functions, respectively. Consequently, the axial u 0 (x 1 ) and transverse w 0 (x 1 ) displacement fields can be written as: where the localized displacement coordinate vectors {U l (x 1 )} and {W l (x 1 )} corresponding to the twonoded 1D element enclosing the point x 1 are given as: and the associated shape function matrices are: The superscripts (·) (1) and (·) (2) in the above equation denote the local node numbers for the two-noded element. Using the definition of the RC derivative in Eq. (4), the fractional derivative D α x1 [u 0 (x 1 )] is expressed as: where D 1 (·) is the first integer-order derivative and s 1 is a dummy variable along the direction x 1 used within the definition of the fractional-order derivative. The above expression can be recast as: is the kernel of the fractional derivative. Note that kernel A(x, y, l, α) is a function of the relative distance between the points x and y, and can be interpreted as an analog form of the attenuation functions typically 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. Similar expressions for can also be derived using an equivalent approach.
). These integer-order derivatives at s 1 are evaluated in terms of the nodal variables corresponding to the element Ω e p , such that s 1 ∈ Ω e p . Using Eq. (19), the integer order derivatives can be expressed as: (23c) Using the above expression for the integer-order derivative D 1 [u 0 (x 1 )], the fractional-order derivative in Eq. (21) can be written as: where the following definition for the attenuation function A p (x 1 , s 1 , l A , l B , α) is employed: Note that the evaluation of the fractional derivative requires a convolution of the integer-order derivatives across the horizon of nonlocality (x 1 − l A , x 1 + l B ). Thus, while obtaining the FE approximation in Eq. (24a), the nonlocal contributions from the different finite elements in the horizon must be properly attributed to the corresponding nodes. In order perform this mapping of the nonlocal contributions from elements in the horizon, we transform the nodal values {U l (s 1 )} and {W l (s 1 )} into the respective global variable vectors {U g } and {W g } with the help of connectivity matrices in the following manner: The connectivity matrices [C u (x 1 , s 1 )] and [C w (x 1 , s 1 )] are defined for the axial and transverse displacement vectors, respectively. These matrices are designed such that they are non-zero only if s 1 lies in the domain (x 1 − l A , x 1 + l B ), the horizon of nonlocality for x 1 . It is evident that these matrices activate the contribution of the nodes enclosing s 1 for the numerical evaluation of the convolution integral given in Eq. (24a) [49]. Using the above formalism, Eq. (24a) can be recast as: and [B θ (s 1 )] have been given in Eq. (23).

B. Nonlinear f-FEM model
In the following, we will use the above formalism to obtain the algebraic equations of motion for the geometrically nonlinear response of the fractional-order nonlocal beam by minimizing the total potential energy of the beam. The first variation of the deformation energy U defined in Eq. (14) is given by: Using the above expression for δU along with the nonlocal axial strain derived in Eq. (11) for the Euler-Bernoulli beam theory and the definitions of the stress resultants given in Eq. (17), we obtain the following expression for the statement of the minimum potential energy δΠ = 0: Note that the variations δu 0 (x 1 ) and δw 0 (x 1 ) are independent of each other. Now, by collecting the terms containing δu 0 (x 1 ) and δw 0 (x 1 ), the above weak form is expressed in the form of the following equations: We recast the above equations in terms of of the displacement field variables by using the expressions for the stress resultants given in Eqs. (9) and (17), and obtain the following: where A 11 = Ebh and D 11 = Ebh 3 /12 are the axial and bending stiffness of the beam, respectively. The above equations in the weak form are converted into algebraic equations using the f-FEM model discussed in §III A. Now, by using Eqs. (27)(28) and enforcing minimization, the above integrals reduce to the following algebraic equations: In the above system of equations, the different stiffness matrices are obtained as: and the axial and transverse nodal force vectors are obtained as: The solution of the system of algebraic equations in Eq. (33) provides the nodal generalized displacement coordinates, which can then be used along with the Euler-Bernoulli relations in Eq. (10) to determine the displacement field at any point of the beam. The effect of the geometric nonlinearity is evident from the expressions of the stiffness matrices given in Eq. (34), where it appears that the stiffness matrices are a function of the deformed system configuration. Note that the linear f-FEM model can be obtained by ignoring the contribution of D α x1 [w 0 (x 1 )] in Eq. (34). Further, similar to classical beam models accounting for geometric nonlinearity, a coupling between the axial and bending responses is obtained through the non-zero [K 12 ] and [K 21 ] matrices.
In this study, we solve the nonlinear algebraic equations using a Newton-Raphson (NR) iterative numerical scheme. Similar to classical nonlinear models, the NR procedure for the fractional-order nonlinear equations also requires the evaluation of the tangent stiffness matrix. The procedure to derive the tangent stiffness matrix [K i T ] corresponding to the i-th iteration of the NR method follows from the classical approach: where {X} is a vector containing all the displacement coordinates of the beam and {R} is the residual vector evaluated at the current iteration. The residual {R} at the i-th iteration step is given by: whereK i S denotes the stiffness matrix of the entire system. In order to obtain the solution of the nonlinear equations, we also adopt a load increment procedure wherein the nonlocal response of the beam at each load step is evaluated using the NR iteration: The iterations at each load level are continued until the residual becomes less than a chosen tolerance.

C. Numerical integration for nonlocal matrices
In this following, we provide the details of the numerical scheme used to integrate both the stiffness and the tangent stiffness matrices. Note that the numerical procedure to evaluate the force vector follows directly from classical FE formulations. Hence, we do not provide a complete account of all the steps. Evaluation of the stiffness matrices for the nonlocal system given in Eq. (34) requires the evaluation of the nonlocal matrices [B ] ( ∈ {u, w, θ}). As seen in Eqs. (27)(28), this involves a convolution of the integer-order derivatives with the fractional-order attenuation function A p (x 1 , s 1 , l A , l B , α) over the nonlocal horizon (x 1 − l A , x 1 + l B ). Clearly, the numerical model for fractional-order derivative involves additional integration over the nonlocal horizon to account for nonlocal interactions. These nonlocal interactions across the horizon of nonlocality have already been accounted for numerically in [18,54]. However, differently from these studies, the attenuation function A p (x 1 , s 1 , l A , l B , α) in the fractionalorder model involves an end-point singularity (more specifically, at x 1 = s 1 ) due to the nature of the kernel (see Eq. (24b)). The fractional-order nonlocal interactions as well as the end-point singularity were addressed in [49] and are briefly reviewed here.
In the following, we describe the procedure for the numerical evaluation of the nonlinear stiffness matrix [K 12 ]. The same procedure can be extended to evaluate the other stiffness matrices, hence their discussion is not provided here. In order to perform the numerical integration of the nonlocal and nonlinear stiffness matrix [K 12 ], we adopt an isoparametric formulation and introduce a natural coordinate system ξ 1 . The Jacobian of the transformation x 1 → ξ 1 is given as J(ξ 1 ). Now, by using the Gauss-Legendre quadrature rule, the matrix [K 12 ] is approximated as: where ξ i,j 1 is the j−th Gauss integration point in the i−th element,ŵ j is the corresponding weight for numerical integration, and N GP is the total number of Gauss points chosen for the numerical integration, such that j ∈ {1, ...N GP }. J i is the Jacobian of the coordinate transformation for the i−th element. The hat symbol on the weight is used to distinguish it from the transverse displacement. As previously highlighted, the matrices [B ] ( ∈ {u, w, θ}) (equivalently, [B (ξ i,j 1 )] in the discretized form) involve an additional integration (see Eq. (27b)) due to the fractional-order nonlocality and are evaluated in the following manner: where x i,j 1 is the Cartesian coordinate of the Gauss point ξ i,j 1 and [B (s 1 )] is given in Eq. (23). Using Eq. (24b), we obtain the following expression for B (x i,j 1 ) : where, the integrands I L and I R are given as: Note that the limits of the integrals in Eq. (40a) span over the elements that constitute the nonlocal horizon ( These limits for the convolution integral follow from the terminals of the RC fractional derivative defined for the study of nonlocal continuum in Eq. (4). These integrals are evaluated numerically in the following manner: Singularity at x i,j are truncated in order to account for asymmetric horizon lengths. This is essential to satisfy frame-invariance as discussed in §II.
As discussed previously, due to the nature of the kernel of the fractional-order derivative, an end-point singularity occurs in the integrals at the Gauss point x i,j 1 in the element Ω e i . This is evident from the definitions of the left and right integrals given in Eq. (41). Following [49], this end-point singularity is circumvented by an analytical evaluation of these integrals over the elements containing the singularities. This analytical evaluation can be carried out by using the expressions for [B ] ( ∈ u, w, θ) given in Eq. (27b). The integrals over the remaining elements (i.e. those without singularities) are evaluated using the Gauss-Legendre quadrature method. The expression for this quadrature based integration corresponding to the nonlocal contribution of the r−th element in the horizon of nonlocality of the Gauss point x i,j 1 is given by: where x r,k 1 is the Cartesian coordinate of the k−th Gauss point in the r−th element along x 1 ,ŵ k is the corresponding weight, and J r is the Jacobian of the transformation for r−th element. The above numerical procedure is schematically illustrated in Fig. 3. We highlight that A p is a function of the Cartesian coordinates, therefore global coordinates should be used in its evaluation. As mentioned previously, the procedure followed for the evaluation of the [K 12 ] stiffness matrix extends directly to the other stiffness matrices defined in Eq. (34), as well as the tangent stiffness matrix defined in Eq. (35). Note that all these matrices are pre-assembled in the global form due to the use of the connectivity matrices in Eq. (26). The resulting matrices are then used in the NR iteration procedure outlined in Eq. (37) to obtain the generalized displacement coordinates for the geometrically nonlinear response of the fractional-order nonlocal beam.

IV. APPLICATION AND NUMERICAL RESULTS
We use the f-FEM developed in §III to study the geometrically nonlinear and fractional-order nonlocal elastic response of isotropic beams. More specifically, the effect of the fractional order α and of the length scales l A , l B on the nonlocal elastic behavior of the nonlinear beam are studied. For all the numerical studies reported in this section, the aspect ratio of the beam is chosen to be L/h = 100 in order to satisfy the slender beam assumption of the Euler-Bernoulli beam theory. The length of the beam is maintained at L = 1 m, the width is kept constant at unity, and the elastic modulus is set at E = 30 GPa. It is assumed here that that the length scales l A and l B for the isotropic beam are equal to a constant l f for points sufficiently within the domain of the beam. These length scales are truncated appropriately for points close to the beam boundaries as discussed in §II. Note that, although the presented f-FEM is capable of simulating the static response of the fractional-order nonlocal beam under both axial and transverse loads, here we present only the results for transverse loads. Both the case of a uniformly distributed load per unit length (UDL) (in N/m) applied along the length of the beam, and of a point load (PtL) (in N) applied at a specified point.
In the following, we first present the results of the validation and convergence study and then we present the static response of the fractional-order nonlocal beam for the various loading conditions. For all the numerical results, the mesh discretization was chosen such that l f /l e = 10, where l e is the length of the discretized element. This choice allows for sufficient number of elements to be included in the horizon of nonlocality at any point hence guaranteeing to accurately capture the nonlocal interactions.  [56], and (b) linear nonlocal continuum (α = 0.9, l f = L/10) compared with [49].

A. Validation
We validated the performance and evaluated the accuracy of the f-FEM by using the following strategies: validation #1: the f-FEM model was solved for α = 1.0, corresponding to the classical elastic solid, and the results were compared with those available in the literature [56]; validation #2: the f-FEM model was solved for non-integer values of α while ignoring the nonlinear effects in Eq. (34), and the results were compared against [49]; and validation #3: the f-FEM results were compared with the exact solution for a nonlinear and fractional-order elastic clamped-clamped beam subjected to a distributed transverse load. The validations #1 and #2 can rather be considered as checks on the robustness of the nonlinear f-FEM. Specific details of the three validation strategies are provided here below.
Validation #1: in this study, the fractional order is taken equal to α = 1.0. This assumption reduces the nonlocal constitutive model (strain-displacement relations) in Eq. (8) to their integer-order counterparts for classical elasticity. The nonlinear f-FEM is solved for this choice of α, l f being arbitrary for the local elastic solid, and the numerical results are compared with [56] in the Fig. (4a). This comparison is performed for beams corresponding to two different boundary conditions, namely, pinned-pinned and clamped-clamped. As evident from Fig. (4a), the match between the transverse displacement obtained numerically from the f-FEM and the results in the literature is excellent. The error between them is less than 1% for all cases.
Validation #2: in this study, the nonlinear effects are ignored for the sake of a comparison with the literature over fractional-order nonlocal and linear elastic response [49]. This simplification is achieved by ignoring the term D α x1 [w 0 (x 1 )] in the definition for stiffness matrices given in Eq. (34), which contributes to the nonlinearity. The resulting linearized f-FEM model was solved for α = 0.9 and l f = L/10 and the numerical results are compared with [49] in the Fig. (4b). Similar to validation #1, this comparison is also performed for two different boundary conditions. As evident from Fig. (4b), the match between the transverse displacement obtained numerically from the f-FEM and the results available in the literature is excellent and the relative error is less than 1% for all cases.
Validation #3: in this study, the following displacement field at the mid-plane of the beam is assumed: Note that the above displacement fields satisfy the clamped-clamped boundary conditions for the beam given in Eq. (16). Although, the above displacement fields are independent of the fractional constitutive parameters α and l f , the axial and the transverse forces required to generate the above elastic response of the beam would be functions of these parameters. These distributed loads are obtained from the strong form of the governing equations given in Eq. (16) using the constitutive relations stated in Eqs. (9) and (17). The closed-form expressions of the distributed axial and transverse loads derived for the exact solutions in Eq. (43) are used in the f-FEM model given in Eq. (33). This strategy is similar to the one presented in [49]. The numerical results obtained from the f-FEM for the mid-plane transverse displacement w 0 (x 1 ) are compared with the exact solution. The results, shown in Fig. (5), are in excellent agreement. The error between the numerically obtained f-FEM solutions and the exact solutions is less than 3% for all cases.

B. Convergence
In the following, we present the results of the sensitivity analysis on the FE mesh size. The convergence of the integer-order FEM with element discretization, referred to as the h−refinement, is well established in the literature [57]. In this study, we noted that in addition to the FE mesh size, the convergence of the f-FEM also depends on the fractional order α, on the relative size of l e , and on the nonlocal length scale l f . It appears that the convergence of the f-FEM depends on the strength of the nonlocal interactions across the nonlocal horizon. Therefore, sufficient number of elements N inf (= l f /l e ) should be available in the influence zone of the fractional-order (nonlocal) interactions. Additionally, mesh refinement across the length of the 1D beam would increase the spatial resolution and hence, result in lower inconsistencies due to the truncation of the nonlocal horizon caused by the ceil and floor operations connected to Eq. (41). This approach would provide a better numerical approximation for the nonlocal matrices, [B ] in Eqs. (40a)- (41). Therefore, convergence of the f-FEM is expected when the number of elements in the influence zone N inf , referred to as the dynamic rate of convergence [49,55], is sufficient to accurately capture the nonlocal interactions. In the context of the current study, we establish the convergence of the f-FEM for two different boundary conditions. The fractional-order parameters are chosen as α = 0.8 and l f = L/10. The results are presented in Fig. (6). In both cases, Fig. (6a) for the clamped-clamped and Fig. (6b) for the pinnedpinned beams, the dynamic rate of convergence is noted to be N inf ≥ 5. Additionally, a comprehensive comparison of the convergence of the f-FEM for various values of the fractional parameters α and l f is presented in Table I. In this table, the convergence may be assessed by comparing the values along a prescribed column on the far right, corresponding to the transverse displacement for a chosen set of fractional constitutive parameters. Table I that for higher values of l f and lower values of α, that is in the case of stronger fractional-order nonlocal interactions, larger number of elements are necessary for the convergence of f-FEM. From the tabulated values, it can be seen that N inf = 10 provides a satisfactory numerical convergence yielding a difference smaller than 1% between successive refinements. Therefore, the f-FEM simulations presented below will be carried out using this mesh discretization.

C. Effect of the fractional-order constitutive parameters
In this section, we use the f-FEM to analyze the static response of the fractional-order nonlocal and geometrically nonlinear Euler-Bernoulli beam. More specifically, we analyze the effect of the fractional parameters α and l f on the response of the beam subject to different loading and boundary conditions.
We start by considering a beam that is clamped at both ends and is subject to a UDL of magnitude q 0 (in N/m). The nonlinear load-displacement curves of the clamped-clamped beam are obtained and compared for different values of α and l f . The results of this comparison are provided in Fig. (7). The effect of the fractional order α, with l f being held constant, is studied in Fig. (7a). Recall that for α = 1.0,  the differential and algebraic form of the governing equations and boundary conditions, reduce to those of a geometrically nonlinear classical (local) elasticity study. Additionally, we also compare the effect of the horizon of nonlocality over the elastic response of the beam in Fig. (7b). As evident from Fig. (7), the transverse displacement increases with increasing degree of nonlocality, that is by reducing the value of α or by increasing the horizon length l f /L. These trends point to a reduction in the stiffness of the fractional-order nonlocal structure. As expected, a convergence of the fractional-order model to the local elastic response is noted for α close to 1.0 and l f /L << 1.  We extend the above study to analyze a beam with both ends pinned. The effect of the fractional model parameters over the nonlocal elastic response of the geometrically nonlinear pinned-pinned beam are shown in Fig. (8). From these figures, it may be noted that the transverse displacement increases while reducing the value of α and increasing the nonlocal horizon length l f . These observations for the pinned-pinned beam agree very well with those noted previously for the clamped-clamped beam. The consistency of these observations regardless of the nature of the boundary conditions indicates the robustness of the fractional order model and marks a net departure from the paradoxical results obtained in the literature for more traditional nonlocal beam theories such as, for example, differential or integrodifferential models of nonlocal elasticity [30,32,58]. This robustness of the fractional-order approach to nonlocality is a direct consequence of the positive definite and self-adjoint nature of the fractional-order formulation [49]. Finally, the nonlocal response of the clamped-clamped and pinned-pinned beams are compared for different loading conditions and the results are shown in Fig. (9). The linear and nonlinear load-displacement curves of the Euler-Bernoulli beam subject to an UDL and with fractional constitutive parameters fixed at α = 0.8 and l f = L/10 are compared in Fig. (9a). A similar comparison for the beam subject to a PtL applied at mid-length is shown in Fig. (9b). It appears that the effect of the geometric nonlinearity is significant for the nonlocal pinned-pinned beams, irrespective of the loading conditions. Note that this is in agreement with the observations from the classical theory of elasticity [56].

V. CONCLUSIONS
The current study laid the foundation for the analysis of geometrically nonlinear and nonlocal solids via fractional-order operators while providing accurate numerical schemes for their solution. More specifically, the geometrically nonlinear formulation for fractional-order nonlocal Euler-Bernoulli beams was developed and its solution was achieved via a dedicated finite-element-based numerical method. The general fractional-order continuum theory was extended to develop a frame-invariant and dimensionally consistent nonlinear fractional-order strain-displacement theoretical framework. This theory is well suited to capture scale effects, nonlocality, and heterogeneity in both complex solids and interfaces undergoing large deformations. The constitutive model for the nonlinear response of fractional-order beams was derived following variational principles. The 1D beam model lays the foundation for a more general nonlinear analysis of fractional-order nonlocal solids. Owing to the complexity of the nonlocal models, a fractional-order FEM (f-FEM) was also developed by means of potential energy minimization. This nonlinear f-FEM provides an efficient and accurate computational approach to obtain the numerical solution of nonlinear integro-differential governing equations with a singular kernel. The robustness of this numerical approach for the solution of nonlocal elasticity problems was demonstrated by showing the consistency of the predicted nonlocal response across different loading and boundary conditions. This latter aspect is in contrast with other traditional approaches to nonlocal elasticity that, occasionally, provide paradoxical results depending on the selection of both loading and boundary conditions. Finally we note that, although the f-FEM was used in this study to solve nonlinear beam models, the same computational framework can be easily extended either to the simulation of systems with higher dimensions or to other nonlinear physical processes described by fractional differ-integral equations.

VI. APPENDIX
In the following, we provide the detailed steps for obtaining the strong form of the governing differential equations given in Eq. (16). The governing equations are obtained by minimizing the potential energy given in Eq. (15). The minimization of the potential energy results in two independent integrals which are given in Eq. (31) because of the independent nature of the variations δu 0 (x 1 ) and δw 0 (x 1 ). It follows from Eq. (31a) that: where x1−l B I 1−α x1 and x1 I 1−α x1+l A are the left and right Riesz integrals [59]. The simplification of Eq. (31) to Eq. (44a) is performed using integration by parts and the definitions for the fractional derivatives. The specific steps leading to the above simplification can be found in [49]. The above discussed simplification extends directly to other integrals in Eq. (31). Combining all the simplifications yields the first variation of the potential energy as: Axial Governing Equation Transverse Governing Equation where D 1 x1 (·) denotes the first integer-order derivative with respect to x 1 . The governing differential equations and the associated boundary conditions given in Eq. (16) are now obtained from the above Eq. (45) by using first principle of variational calculus. Note that, for longitudinal boundaries of the beam, i.e., x 1 = {0, L}, l A and l B are equal to zero in the Caputo and RL derivatives defined in Eqs. (4,18), respectively, and the Riesz integrals given in Eq. (44b). These limit cases reduces the boundary conditions to the integer-order counterparts as given in Eq. (16) [40,49].

VII. ACKNOWLEDGEMENTS
The following work was supported by the Defense Advanced Research Project Agency (DARPA) under the grant #D19AP00052, and the National Science Foundation (NSF) under the grant DCSD #1825837. 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.