Non-standard and constitutive boundary conditions in nonlocal strain gradient elasticity

Zaera et al. (Int J Eng Sci 138:65–81, 2019) recently showed that the nonlocal strain gradient theory (NSGT) is not consistent when it is applied to finite solids, since all boundary conditions associated to the corresponding problems cannot be simultaneously satisfied. Given the large number of works using the NSGT being currently published in the field of generalized continuum mechanics, it is pertinent to evince the shortcomings of the application of this theory. Some authors solved the problem omitting the constitutive boundary conditions. In the current paper we show that, in this case, the equilibrium fields are not compatible with the constitutive equation of the material. Other authors solved it omitting the non-standard boundary conditions. Here we show that, in this case, the solution does not fulfil conservation of energy. In conclusion, the inconsistency of the NSGT is corroborated, and its application must be prevented in the analysis of the mechanical behaviour of nanostructures.


Introduction
Although the origin of higher order continuum theories can be found in the second part of nineteenth century and in the beginning of twentieth century, the explosive growth of the application of nanostructures in several scientific and technological fields has renewed the interest on developing this kind of approaches. The reasons can be found in: (i) the impossibility to address with the classical approaches (scale-free) the size effect, which commonly appears in nanotechnology applications; (ii) the lower computational cost in comparison to molecular dynamic techniques.
In the last 15 years, the most widely used generalized continuum approaches to address mechanical problems at nanoscale are the modified strain gradient elasticity theory and the nonlocal elasticity frameworks. The modified strain gradient theory [1] is based in previous formulations by Mindlin [2] and Fleck and Hutchinson [3], and only three size-dependent constants are needed for the case of linear isotropic elastic materials. Thai et al. [4] recently published relevant examples of the application of this theory to the analysis of the mechanical behaviour of nanostructures. Moreover, it is worth to cite the works by Morassi and coworkers [5][6][7][8][9] to address problems related with nanosensors. Other widely used approaches to address the mechanical behaviour of nanomechanical systems fall into the nonlocal continuum mechanics framework. The bases of the theory were formulated by different authors in the 60s of the twentieth century [10][11][12]. Some time later, an integral constitutive model valid for linear isotropic elastic materials was proposed by Eringen et al. [13][14][15][16]. The main characteristics of the so called strain-driven nonlocal elasticity is that the stress at any point of the domain depends on the strain at each point in the solid through a convolution integral with a smoothing kernel. The nonlocal constitutive equation can be written in differential form, if a Helmholtz-type kernel is used in the integral formulation [16]. Based on this transformation, the differential approach has been widely used in a very large numbers of papers, starting with the work of Peddieson et al. [17]. The reader can found more details on the application of this differential model in the reviews recently published by Eltaher et al. [18], Rafii-Tabar et al. [19], and Thai et al. [4]. However, Romano et al. [20] recently pointed out that the pure nonlocal strain-driven elasticity theory leads to problems with no solution in general, due to the ill-posedness of the model. The reason is the requirement to accomplish additional boundary conditions related to the integral nature of the constitutive equations. To overcome these drawbacks, several alternatives have been developed inside the framework of nonlocal theories. Romano et al. [21] proposed the so-called stress-driven nonlocal model, in which the elastic strain at a certain point in the solid is related to the stress at all points of the domain by a convolution integral with a smoothing kernel of the Helmholtz-type. This approach allows to obtain consistent solutions corresponding to the mechanical behaviour of several kinds of nanostructures [22][23][24][25][26]. The mentioned ill-posedness affecting the pure straindriven nonlocal approach can be also avoided using the two-phase local/nonlocal strain-driven constitutive approach, proposed by Eringen [13,27] and recently applied by different authors to analyse both the static and dynamic behaviour of nanostructures [28][29][30][31][32]. Moreover, Barretta et al. [33,34] have recently developed the two-phase local/nonlocal stress-driven formulation.
Lim et al. [35] combined both the fully nonlocal elasticity theory of Eringen and the strain gradient elasticity in a single framework called nonlocal strain gradient elasticity theory (NSGT). The new formulation allows to analyse stiffer (or more flexible) structures than those corresponding to the classical case, appropriately selecting the two scale-dependent material parameters included in the model. The theory has become very popular, and nowadays it is being currently used by a large number of authors, see [36][37][38][39][40][41][42][43][44][45] just to mention a few. When it is applied to bounded domains, the fulfilment of both standard and non-standard boundary conditions, together with the constitutive boundary conditions are required. In some of the quoted papers the solution only satisfies the standard and the non-standard boundary conditions, dismissing the constitutive ones [36,37,38,39,42,45], meanwhile in other works only satisfies the standard and the constitutive boundary conditions, dismissing the non-standard ones [40,41,43,44]. In this respect Zaera et al. [46] showed the inconsistency of the NSGT when applied to finite solids. The requirement to meet the three types of boundary conditions at once, standard, non-standard, and constitutive, makes the problem overconstrained, with no solution in general. Yet, given the large number of works using the NSGT being currently published in the field of generalized continuum mechanics, it is pertinent to evince the shortcomings of the application of this theory. Indeed, in a recent paper [44] it is claimed that, if the NSGT is properly formulated in variational terms, only the standard and the constitutive boundary conditions are required, thus amending the quoted inconsistency of the theory.
In order to stress the aforementioned inconsistency, in the present work we point out the consequences of dismissing certain mandatory boundary conditions, either non-standard or constitutive, to derive an apparently well-posed problem from a factually overconstrained problem. In particular, we analyse the bending behaviour of a simply supported Bernoulli-Euler beam subjected to an uniformly distributed static load, which was considered in different works [39,44,46]. Through this straightforward example we show that: (i) if constitutive boundary conditions are omitted, the equilibrium fields are incompatible with the constitutive equations of the material, as it was pointed out by Barretta and de Sciarra [40]; (ii) if non-standard boundary conditions are omitted, the energy conservation principle is infringed. Therefore the NSGT is inconsistent because, in general, there is no solution if it is applied to bounded solids. Consequently, its use must be definitely prevented to assess scale effects in mechanical nanosystems.
The paper is organised as follows. Section 2 presents the formulation of the studied problem. In Sect. 3, the consequences of dismissing the constitutive boundary conditions are highlighted, whereas Sect. 4 points out those of disregarding the nonstandard ones. Finally, Sect. 5 stresses the main conclusion of this work.

Static bending of a NSGT Bernoulli-Euler beam
The essential hypothesis for the constitutive behaviour of the NSGT considers that the internal elastic energy is composed by two terms: one accounting for the nonlocal effects of the classical strain field and another one accounting for the nonlocality of the strain gradient field [35]. From this assumption and applying a classical variational principle it is possible to identify the total stress as the difference between the classical stress (work conjugate of classical strain) and the gradient of the higher order stress (work conjugate of strain gradient).

Governing equation, standard and nonstandard boundary conditions
Let us consider a simply supported beam of length L and uniform section area A, inertia I, and Young modulus E, subject to a distributed transverse load q z . Axial, out-of-plane, and transverse coordinates are represented respectively as x, y, z, and the corresponding displacements as U x , U y , U z . According to the Bernoulli-Euler beam kinematics, we have where w represents the displacement of the section's centroid in the transverse direction and, for simplicity, Á ð Þ 0 denotes the derivative with respect to the spatial variable x. The axial strain, e, is given by Specialising the general NSGT constitutive equation for the one dimensional problem, the normal axial stress is given by [35,46] where r 0 and r 1 represent nonlocal and higher order stresses, respectively defined as ð4Þ l and j being parameters accounting for scale effects related to the strain gradient field e 0 and to the nonlocal stress field respectively. The kernel appearing in above equations is given by Assuming the previous hypotheses, the bending moment is given by with M 0 ¼ ÀEI and being respectively nonlocal and higher order moments. According to Lim et al. [35], and particularizing the general 3D formulation to the 1D case, the internal energy density potential (per unit volume) U is given by the product of work conjugate pairs r 0 and e, on one side, and r 1 and e 0 , on the other From the kinematics of the Bernoulli-Euler beam model, the internal energy of the whole beam is derived as with Due to the conservative character of the system, governing equation and boundary conditions are derived from the principle of minimum total potential energy dP ¼ dU À dV ¼ 0, where P and V are the total potential energy and the potential of applied forces, respectively. This principle implicitly imposes the equivalence between the work of external forces and the internal energy of the solid, implying conservation of energy. Now, integrating by parts and considering the symmetry of the 1D kernel [Eq. (6)] with respect to axial coordinates x and x, we get the first variation The first variation of dV is given by By virtue of the principle of minimum total potential energy [47], we get the equilibrium equation as well as the standard essential and natural boundary conditions to be imposed at as well as non-standard ones It should be noted that non-standard boundary conditions (20)-whether essential or natural-have to be fulfilled, together with the standard boundary conditions (18) and (19), to achieve equilibrium (dP ¼ 0). Now introducing the following nondimensional variables and parameters the dimensionless governing equation is given by Standard and non-standard boundary conditions to be imposed at n ¼ 0 and n ¼ 1 are written as w ¼ 0; or where now Á ð Þ 0 represents derivative with respect to the nondimensional spatial variable n. For the study case (simply supported beam), the standard boundary conditions read w 0 Non-standard essential boundary conditions, i.e. w 00 0 ð Þ ¼ w 00 1 ð Þ ¼ 0, have been selected as in [39]. The integral constitutive equations read where nonlocal moment is given by and higher order moment by

Constitutive boundary conditions
The equivalence between the integral and differential constitutive equations, the last one equipped with constitutive boundary conditions, can be set as follows [48]. The linear Fredholm integral equation of first kind with exponential kernel is equivalent to a boundary value problem consisting on the ordinary differential equation As stated by Romano et al. [20], the quoted integral equation implicitly contains the boundary conditions which are necessary and sufficient for the existence and uniqueness of its solution.
The constitutive equation (27) relates the nonlocal moment with the curvature through the integral equation (29), with a ¼ 0, b ¼ 1, l ¼ À1=h, u ¼ w 00 s ð Þ, and f ¼ À2h M 0 , thus being equivalent to the differential constitutive equation subject to the constitutive boundary conditions Likewise, the integral constitutive equation (28) relates the higher order moment with the derivative of the curvature through the integral equation (29), Þ, and f ¼ À2h=g 2 M 1 , thus being equivalent to the differential constitutive equation subject to the constitutive boundary conditions The differential equation relating the total moment and the curvature is derived from the two previous BVPs. From Eqs. (26), (32), and (34) we get the constitutive equation Differentiating Eq. (26) and considering Eq. (34) Particularizing Eq. (37) for n ¼ 0 and using BCs (33) and (35), we get Similarly for n ¼ 1, we get The constitutive boundary conditions (38) and (39) are necessary and sufficient for the existence and uniqueness of the solution of the integral constitutive relation (26) [40]. Similarly, the integral constitutive equation (26) is equivalent to the differential constitutive equation (36) equipped with constitutive boundary conditions (38) and (39).

Formulation in the displacement variable
Now the problem will be formulated in terms of the displacement w. From Eqs. (36) and (22), the total bending moment follows the expression M ¼ Àh 2 q z À w 00 þ g 2 w IV : Taking now M 00 from Eq. (40), the governing Eq. (22) can be written as The previous equation is subject to the following boundary conditions: • Four standard boundary conditions: For the reasons stated above, the previous eight boundary conditions are all of them mandatory, and the sixth-order governing Eq. (41) has more boundary conditions than those required. In the following sections, we outline the serious implications of dismissing the constitutive boundary conditions, as well as those of not using the non-standard ones. To that aim, we consider the bending behaviour of a simply supported Bernoulli-Euler beam subjected to an uniformly distributed static load q z ¼ q, which was considered in different works [39,44,46] as pointed out in the introduction section. After that, the inconsistence of the NSGT clearly emerges.

Consequences of dismissing the constitutive boundary conditions
The general solution of Eq. (41) for uniformly distributed load ( q z ¼ q) is given by A 0 to A 5 being six arbitrary constants that have to be determined imposing six boundary conditions. First, the problem is solved with the four classical boundary conditions ( , and two non-standard boundary conditions ( w 00 0 ð Þ ¼ w 00 1 ð Þ ¼ 0), thus disregarding the two constitutive boundary conditions (42) and (43). In this case, and after some algebra, the displacement field is given by which provides the same results presented in Table 2 of [39] for the parameters considered in that work. The displacement field does not fulfil the constitutive boundary conditions since both left terms of expressions (42) and (43) are non-nil in general qðg À hÞðg þ hÞ 2g tanh 1 2g À 1 Thus the bending moment field derived from the constitutive equation is different to that derived from the internal equilibrium equation, as it will be shown next.
Using expressions (27) and (28), the nonlocal and higher order bending moments read, respectively Finally, the bending moment field derived from the nonlocal strain gradient constitutive equation is given by On the other side, due to the isostatic character of the considered beam, the distribution of bending moments can be also determined by integration of the internal equilibrium Eq. (22) subject to the standard boundary conditions M 0 ð Þ ¼ 0 and M 1 ð Þ¼ 0, leading to Therefore, the bending moment field derived from the constitutive equation, given by expression (49), is incompatible with that derived from the equilibrium equation. The mismatch is also evident in Fig. 1, and is the consequence of dismissing the constitutive boundary conditions in the integration of the governing equation.

Consequences of dismissing the non-standard boundary conditions
Now, the sixth-order governing Eq. (41) is solved with the four classical boundary conditions ( , and the two constitutive boundary conditions (42) and (43) as presented in both [46] and [44], thus disregarding the two non-standard boundary conditions. It has to be noted that the solutions derived herein are fully equivalent to that considered in [44] when the mixture parameter a used in the quoted work is set to zero to recover the pure NSGT theory. Nevertheless, the inconsistency that will apear for a ¼ 0 still holds for a mixed nonlocal strain gradient formulation with a 6 ¼ 0.
In the present case (a ¼ 0), the displacement field is given by which does not satisfy the non-standard boundary conditions, neither essential nor natural, as shown in [46]. Equation (51) provides results coincident to those presented in Tables 3 and 4 of [44] for the parameters considered in that work. Using expressions (27) and (28), the nonlocal and higher order bending moments read, respectively It can be shown that now the bending moment field M ¼ M 0 À M 0 1 is consistent with that derived from the internal equilibrium equations, given by Eq. (50). However, a deep inconsistency still arises when nonstandard boundary conditions are dismissed, as it will be discussed next.
External work in the beam is given by On the other side, the internal energy is calculated as (see Eqs. (11)-(13) in this paper, or Eq. (10) in [44]) U ¼ 1 2 where the minus sign in the integrands of Eq. (55) appears due to the specific sign criterion considered here. The previous expressions (54) and (55) show that, in general, W ext and U are not equal (see Fig. 2 for a comparison of both expressions for a particular case with q ¼ 1). The implications of the mismatch between external work and internal energy are Nondimensional Energies (c) Fig. 2 Variation of the external work W ext , of the strain energy U, and of the work due to the boundary term W B with the nonlocal parameter h, for a simply supported beam with uniformly distributed load q ¼ 1. a g ! 0; b g ¼ 0:1; c g ¼ 0:2. Solution derived dismissing non-standard boundary conditions. All energies have been nondimensionalized with the external work corresponding to the classical (local) case, profound, since the conservation of energy is infringed by a solution that only fulfils the standard and constitutive boundary conditions, dismissing nonstandard ones.
Recently, in the paper by Barretta and Marotti de Sciarra [44] it is claimed that no ill-posedness holds when dismissing non-standard boundary conditions if the NSGT is established by an adequate variational formulation with appropriate test fields. This approach is posed through the application of the principle of virtual work (Eq. (1) in [44]), being fully equivalent to the principle of minimum total potential energy for elastic solids. This principle is stated through the definition of the internal work as where here we are considering the work variables integrated along the loading process. This definition of the internal work, valid for the classical bending problem, is not applicable to NSGT beams. The correct one has to be derived as the product of conjugate pairs as in Eq. (11) [or Eq. (55)]. Therefore, in order to sort out this shortcoming, these authors impose a second variational condition (Eq. (11) in [44]) stating the equivalence between internal work, defined as in (56), and internal energy, defined as in Eq. (55) [or in Eq. (10) of [44], taking the mixture parameter a ¼ 0]. When the first variation of the internal energy is further developed, a boundary term appears [see Eq. (12) of that paper]. This term is then dismissed by choosing a virtual field of curvatures taking null value at the boundaries, d w 00 0 ð Þ ¼ 0 and d w 00 1 ð Þ ¼ 0. However, it has to be pointed out that selecting a virtual field with null values at n ¼ 0 and n ¼ 1 implies that the primary variable w 00 has to be specified at these points, i.e. essential boundary conditions are imposed [49]. Alternatively, the boundary term could be cancelled by imposing the natural boundary conditions M 1 0 ð Þ ¼ 0 and M 1 1 ð Þ ¼ 0, or any other suitable choice of essential and natural boundary conditions. All in all, from the nature of the variational problem, the fulfillment of the surplus non-standard boundary conditions, not considered in the quoted paper, are required. The verification of this statement comes from the integration along the loading process of the boundary term The addition of this term to the external work now equals the strain potential, i.e., W ext þ W B ¼ U, as it can be deduced from Eqs. (54), (55) and (57) (see Fig. 2).
There are specific cases where the equivalence W ext U holds. The first one, of limited conceptual interest for the current discussion, is that corresponding to h ¼ g, for which both energies are coincident with those derived for the classical (local) beam (Fig. 2). In this case, the stiffening effect of the gradient part appears to offset the softening effect of the nonlocal part, and the displacement field equals that of a classical supported beam. The second case is of much greater interest and deserves to be commented. If the nonlocal parameter h tends to zero, the energy balance is also achieved for any value of the gradient parameter g. Here the nonlocal kernel becomes the Dirac delta and both the constitutive and the governing equations drop to the strain gradient ones, for which non-standard boundary conditions are also required [50]. It may appear contradictory that now a solution obtained without imposing mandatory non-standard boundary conditions does satisfy the energy balance. The rationale is found in the imposed constitutive boundary conditions (42) and (43) leading to w 000 0 ð Þ ¼ 0 and w 000 1 ð Þ ¼ 0 when h ! 0, which are coincident with the non-standard natural boundary contitions M 1 ð0Þ ¼ 0 and M 1 ð1Þ ¼ 0 (notice that, if the nonlocal parameter tends to zero, M 1 ¼ Àg 2 w 000 , see Eq. (82) in [46]). As a particular case, the classical (local) formulation recovered for h ! 0 and g ! 0 also satisfies the energy balance between external work and internal energy.

Conclusion
The mandatory character of the standard boundary conditions, as well as of the non-standard and constitutive ones, leads to an overconstained problem with no solution in general, as stated in [46]. In this work it is showed that, if constitutive boundary conditions are dismissed, the equilibrium fields are not compatible with the constitutive equations of the material. Moreover, we pointed out that if nonstandard boundary conditions are omitted, the energy conservation principle is infringed. Therefore, the inconsistency of the NSGT when applied to finite structures is corroborated, and its use to assess scale effects in nanostructures must be definitely prevented.
However, there are other fully consistent nonlocal models: the two-phase local/nonlocal strain driven elasticity [13,29,32], the nonlocal stress-driven elasticity [21], or the two-phase local/nonlocal stress-driven elasticity [33]. All of them emerge as alternatives to the NSGT.

Remarks
A typo detected in the expression of the displacement field w in Section 6.3 of [46] has been corrected in Eq. (51) of the current manuscript. However, the expression of the curvature and the values presented in Fig. 3a of the quoted work are correct.
All the closed-form expression in the present communication have been obtained with the Mathematica programming environment [51].