Nonlocal Mechanical Behavior of Layered Nanobeams

Abstract

The research at hand deals with the mechanical behavior of beam-like nanostructures. Nanobeams are assembled of multiple layers of different materials and geometry giving a layered nanobeam. To properly address experimentally noticed size effects in structures of this type, an adequate nonlocal elasticity formulation must be applied. The present model relies on the stress-driven integral methodology that effectively circumvents known deficiencies of other approaches. As a main contribution, a set of differential equations and boundary conditions governing the underlaying mechanics is proposed and applied to two benchmark examples. The obtained results show the expected stiffening nonlocal behavior exhibiting most of smaller and smaller structures and modern devices.

1. Introduction

The mechanical behavior at the micro- and nanoscale is governed by different forces and potentials than the macroscale structures [1]. It is well-established that relations between forces and displacements are influenced by small dimensions of such structures, nowadays widely employed to model nanocomposites [2,3] and new-generation devices, such as micro-bridges [4], nano-switches [5], nano-generators [6], and energy harvesters [7,8,9]. In other words, the whole neighborhood of a point contributes to these relations and the term nonlocal is appropriately used in these circumstances. For instance, early experimental observations [10] find that silicon is strongly influenced by size-effects. Measurements show that Young’s modulus varies from approximately 80 GPa for lower nanocantilever thicknesses, converging to 170 GPa for higher thicknesses. The upper limit value corresponds to $\left[110\right]$ bulk silicon. For the investigation at hand, beams assembled of several layers of different materials are of particular interest. Thin metal films are typically used in such applications. Bending load is often introduced by means of contactless electrostatic attraction. In [11], an experimental analysis of a nanocantilever made of copper and silicon layers is presented. It is found that the rather high Young’s modulus of copper (approximately 750 GPa) for low thicknesses decreases dramatically with the increase in thickness, converging to about 120 GPa. Hence, to analyze structures of such small dimensions, one needs a different mechanical apparatus than the macroscale engineer.

The aforesaid experimental observations triggered the high interest of researchers of different specialties. A detailed insight into the size-effects of nanostructures can be achieved by the molecular dynamics (MD) [12] or molecular structural mechanics [13,14,15]. Such simulations can be computationally intensive and, at least for applications involving sensors, analytical nonlocal beam models are preferred. To complicate the issue even further, most of the results dealing with the laminated beams fall within the scope of dynamics and cover a variety of boundary conditions, geometry, and loadings. As an illustration, several examples are mentioned. In [16], a nonlinear finite element analysis of the bending of a rotating laminated nanocantilever is analyzed. The vibrational framework for electro-magneto-elastic bending analysis of three-layer nanobeams is presented in [17]. The differential formulation involving the stress gradient is used. The formulation is extended to curved nanobeams in [18]. Attempts to include temperature effects can also be found (see [19,20]). Feng et al. [20] clearly demonstrated that both the cantilever dimensions and even small temperature changes contribute to the statics and dynamics of layered nanocantilevers.

The influence of the interface and surface effects in layered nanobeams is analyzed in [21] by means of analytical models. Although the approach relies on a different size effect methodology than this paper, some results have to be pointed out. In particular, the interface and surface effects can have in some configurations (especially in beams with a large length/thickness ratio) a profound influence on the elastic deformation; thus, a nonlinear analysis incorporating these effects could be more suitable. An interested reader is referred to the work of Müller and Saúl [22] for a comprehensive analysis on the subject and that of Miller and Shenoy [23] for another approach to size-dependent behavior of nanosized structural elements, nanoplates and nanobeams in particular. When dynamics of coated layered beams of microscopic thickness is to be addressed, then a surface model including damping can also have an important role. For instance, Rongong et al. [24] presented an experimental work about a novel damping layered coating, reporting important properties of the damping layer. Other sources [25,26] report that experimental verification in some cases confirms the importance of the interface effects, while in other cases the interface effects can be disregarded. Thus, one should carefully study the layered beam in question to decide if the surface and interface effects are important or not.

Most of the previous efforts are related to the gradient-based nonlocal formulation. The original theory, as proposed by Eringen [27], is based on the strain-driven integral convolution leading to an equivalent gradient formulation. These results became a popular starting point for many nonlocal applications in mechanics. However, the theory involves some assumptions that are frequently overlooked. In particular, the original contribution was developed to describe surface waves in unbounded continua, thus applications to bounded continua such as beams or plates can be affected by a non-physical behavior. To assure equivalence of the integral and gradient formulation, a suitable set of boundary conditions must be provided. Unfortunately, these boundary conditions contradict equilibrium constraints, as shown in [28,29]. The paradoxical behavior has been known to occur [30,31] for quite some time. This is particularly true for the nanocantilever with a tip loading, a frequently used nanostructure. An efficient answer to this problem is the stress-driven integral formulation [29,32] in which the strain state at a point is influenced by the stress distribution in the neighborhood. It also introduces the constitutive boundary conditions necessary for the correct solution of the equivalent differential formulation. In strong contrast to existing formulations, the present research takes the stress-drive integral formulation as the starting point.

This paper aims to contribute to the current body of literature on the mechanics of small-scale structures in nonlocal integral elasticity. To this end, the stress-driven integral formulation [32] is extended in the present paper to effectively capture scale phenomena in composite beams made of an arbitrary number of layers of different constitutive materials, with an arbitrary width of each layer. The presented model can accommodate beams simultaneously loaded by axial and bending loads.

2. Fundamental Concepts

2.1. Geometry of the Beam

The present section defines assumptions regarding the geometry for beams assembled of n layers of different materials (Figure 1). Bending in the $x-z$ plane is analyzed by means of the Bernoulli–Euler formulation. We select the x-axis as the longitudinal axis, while the vertical position of the origin of the coordinate system is situated at the geometric centroid of the cross-section. Therefore, some layer $i\in \{1,2,\dots ,n\}\subset \mathbb{N}$ is represented by the rectangular cross-section of height $h_i=z_i-z_{i-1}$. It is assumed that coordinates $z_i$ and $z_{i-1}$ define the upper and lower side of layer i. The width $b_i$ is defined in the y direction, i.e., orthogonal to the bending plane. Indices throughout the paper denote the layer number and are not used in the sense of the Einstein summation convention. The beam’s length is denoted by L.

Thus, the cross-section of the beam is assumed to be composed of n rectangles. If the cross-section is symmetric with respect to the $x-z$ plane and if the origin of the coordinate system O is situated at the centroid, then the first moment of area must vanish:

$${\int }_{\mathsf{\Omega }}\overline{z}\mathrm{d}A=\sum _{i=1}^n{b}_i{\int }_{z_{i-1}}^{z_i}\overline{z}\mathrm{d}\overline{z}=\frac{1}{2}\sum _{i=1}^n{b}_i(z_i^2-z_{i-1}^2)=0.$$

(1)

However, due to differences in the Young’s modulus of each layer, the physical neutral surface in which the normal stresses vanish does not correspond in general to the geometrical middle surface [33]. A similar situation occurs in functionally graded beams along the vertical coordinate z (see [34,35,36]). The distance between these two surfaces is denoted by ${\zeta }_0$.

2.2. Notation

To simplify elaborations, new notation for vectors describing geometry and material of layers can be introduced in a more compact manner as:

$$\begin{array}{c}\mathbf{b}=\left\{\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right\},\mathbf{h}=\left\{\begin{array}{c}{h}_1\\ h_2\\ ⋮\\ h_n\end{array}\right\},\mathbf{A}=\mathbf{b}\odot \mathbf{h}=\left\{\begin{array}{c}{b}_1{h}_1\\ b_2{h}_2\\ ⋮\\ b_n{h}_n\end{array}\right\},\\ \mathbf{S}=\left\{\begin{array}{c}{b}_1{\int }_{z_0}^{z_1}(\overline{z}-{\zeta }_0)\mathrm{d}\overline{z}\\ b_2{\int }_{z_1}^{z_2}(\overline{z}-{\zeta }_0)\mathrm{d}\overline{z}\\ ⋮\\ b_n{\int }_{z_{n-1}}^{z_n}(\overline{z}-{\zeta }_0)\mathrm{d}\overline{z}\end{array}\right\},\mathbf{I}=\left\{\begin{array}{c}{b}_1{\int }_{z_0}^{z_1}{(\overline{z}-{\zeta }_0)}^2\mathrm{d}\overline{z}\\ b_2{\int }_{z_1}^{z_2}{(\overline{z}-{\zeta }_0)}^2\mathrm{d}\overline{z}\\ ⋮\\ b_n{\int }_{z_{n-1}}^{z_n}{(\overline{z}-{\zeta }_0)}^2\mathrm{d}\overline{z}\end{array}\right\},\\ {\mathbf{L}}_{\lambda }=L\mathbf{\lambda }=L\left\{\begin{array}{c}{\lambda }_1\\ {\lambda }_2\\ ⋮\\ {\lambda }_n\end{array}\right\}=\left\{\begin{array}{c}{L}_{\lambda ,1}\\ L_{\lambda ,2}\\ ⋮\\ L_{\lambda ,n}\end{array}\right\},\mathbf{E}=\left\{\begin{array}{c}{E}_1\\ E_2\\ ⋮\\ E_n\end{array}\right\},\mathbf{1}=\left\{\begin{array}{c}1\\ 1\\ ⋮\\ 1\end{array}\right\},\end{array}$$

(2)

where the last vector consists of n ones and the symbol ⊙ represents the Hadamard product, i.e., the component-wise multiplication of two vectors. Elements of the vector $\mathbf{b}$ are widths and $\mathbf{h}$ heights of each layer; and $\mathbf{A}$, $\mathbf{S}$, and $\mathbf{I}$ contain cross-section areas and first and second moments of cross-section area for each layer, respectively. Young’s moduli of layers $E_i$ are given as elements of $\mathbf{E}$. The nonlocal parameters that are described below are given in ${\mathbf{L}}_{\lambda }$. Although using boldface uppercase Latin letters violates standard vector notation, its meaning is hopefully more easily interpreted by readers due to familiarity to the equivalent scalar notation.

The following scalar quantities describing local stiffnesses are also introduced

$$\begin{array}{c}{k}_{EA}=\sum _{i=1}^n{b}_i{E}_i{\int }_{z_{i-1}}^{z_i}\mathrm{d}\overline{z}=\sum _{i=1}^n{b}_i{E}_i{h}_i=\sum _{i=1}^n{A}_i{E}_i=(\mathbf{E}\odot \mathbf{A})·\mathbf{1},\hfill \\ k_{ES}=\sum _{i=1}^n{b}_i{E}_i{\int }_{z_{i-1}}^{z_i}(\overline{z}-{\zeta }_0)\mathrm{d}\overline{z}=(\mathbf{E}\odot \mathbf{S})·\mathbf{1},\hfill \\ k_{EI}=\sum _{i=1}^n{b}_i{E}_i{\int }_{z_{i-1}}^{z_i}{(\overline{z}-{\zeta }_0)}^2\mathrm{d}\overline{z}=(\mathbf{E}\odot \mathbf{I})·\mathbf{1},\hfill \end{array}$$

(3)

while the nonlocal contributions are denoted by

$$\begin{array}{c}{k}_{EA}^{NL}=\sum _{i=1}^n{b}_i{E}_i{h}_i{L}_{\lambda ,i}^2=(\mathbf{A}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1},\hfill \\ k_{ES}^{NL}=\sum _{i=1}^n{b}_i{E}_i{L}_{\lambda ,i}^2{\int }_{z_{i-1}}^{z_i}(\overline{z}-{\zeta }_0)\mathrm{d}\overline{z}=(\mathbf{S}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1},\hfill \\ k_{EI}^{NL}=\sum _{i=1}^n{b}_i{E}_i{L}_{\lambda ,i}^2{\int }_{z_{i-1}}^{z_i}{(\overline{z}-{\zeta }_0)}^2\mathrm{d}\overline{z}=(\mathbf{I}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1}.\hfill \end{array}$$

(4)

The symbol · denotes the scalar product of two vectors as usual.

2.3. Kinematics of the Beam

The presence of the shift of neutral surface ${\zeta }_0$ requires the modification of the standard equation representing the longitudinal displacement $u(x,z)$ of each point in the $x-z$ plane:

$$u(x,z)=u_0\left(x\right)+\phi \left(x\right)(z-{\zeta }_0).$$

(5)

The first term is an average longitudinal displacement due to axial loading. Differences in the stiffness of each layer are causing the variation of the longitudinal displacement in the transverse direction even in the simplest case of axial loading. In that way, it is more convenient to work with an average longitudinal displacement. The average longitudinal displacement in a particular cross-section x can be evaluated as:

$$u_0\left(x\right)=\frac{1}{A}{\int }_{\mathsf{\Omega }}u(x,z)\mathrm{d}A=\frac{1}{A}\sum _{i=1}^n{b}_i{\int }_{z_{i-1}}^{z_i}u(x,\overline{z})\mathrm{d}\overline{z},$$

(6)

where $A={\sum }_{i=1}^n{A}_i=\mathbf{A}·\mathbf{1}$ is the total cross-section area. The second term in Equation (5) accounts for the rotation of the cross-section due to bending loads.

The Bernoulli–Euler hypothesis imposes the absence of the shear stresses and thus relates the transversal displacement $w\left(x\right)$ in the direction of the z-axis and angle of rotation $\phi$ in a usual manner:

$$w^{\left(1\right)}\left(x\right)=-\phi \left(x\right),$$

(7)

where the apex ${}^{\left(n\right)}$ denotes the nth derivative with respect to x. Differentiation of the longitudinal displacement with respect to x now gives the normal strain:

$$\epsilon (x,z)=u_0^{\left(1\right)}-w^{\left(2\right)}(z-{\zeta }_0).$$

(8)

3. Nonlocal Material Model

Application of the stress-driven nonlocality within the beam framework was initially introduced in a series of papers [29,32] to overcome applicative difficulties of the strain-driven purely nonlocal approach [28]. The same concept can be adopted for layered beams as well. The normal strain in some layer i depends on the axial stress $\sigma (x,z)$ in the neighborhood, and the closest points have more influence than the more distant ones:

$$\epsilon (x,z)={\int }_0^L{\varphi }_{\lambda ,i}(x-\xi )E_i^{-1}\sigma (\xi ,z)\mathrm{d}\xi .$$

(9)

The kernel function ${\varphi }_{\lambda ,i}(x-\xi )$ can be selected in several ways. In this particular instance, the exponential form is conveniently chosen. The kernel is a function of a single variable—the longitudinal coordinate x—thus only normal stresses in points with the same coordinate z contribute to the nonlocality:

$${\varphi }_{\lambda ,i}\left(x\right)=\frac{1}{2L_{\lambda ,i}}exp\left(-\frac{\left|x\right|}{L_{\lambda ,i}}\right).$$

(10)

The parameter governing nonlocal behavior is the characteristic length $L_{\lambda ,i}={\lambda }_{i}L$, where ${\lambda }_i$ is the small-size parameter of a particular material layer. The reader should be warned that the value of this parameter for a particular material is still unknown in most cases. The parameter should be ideally determined experimentally, but this is rarely done so (see [37,38,39] for some attempts). Alternatively, molecular dynamics or molecular structural mechanics can be employed [15,40,41,42,43,44] for the purpose, but the reported results also indicate that the precise value of the small-size parameter is far from being accurately determined.

The equivalent differential problem can be stated in each layer as [45,46]:

$$\begin{array}{c}\sigma =E_i\left(-L_{\lambda ,i}^2{\epsilon }^{\left(2\right)}+\epsilon \right),\hfill \end{array}$$

(11)

with the following boundary conditions:

$$\begin{array}{c}{\left.L_{\lambda ,i}{\epsilon }^{\left(1\right)}-\epsilon =0\right|}_{\mathrm{at}(0,z)},\hfill \\ {\left.L_{\lambda ,i}{\epsilon }^{\left(1\right)}+\epsilon =0\right|}_{\mathrm{at}(L,z)}.\hfill \end{array}$$

(12)

In the context of beams, the nonlocal framework involving stresses and strains is not that useful. The usual practice is to work with displacements and stress resultants. The normal force in the cross-section follows from the balance of the linear momentum:

$$N\left(x\right)=\sum _{i=1}^n{b}_i{\int }_{z_{i-1}}^{z_i}\sigma (x,\overline{z})\mathrm{d}\overline{z}=\sum _{i=1}^n{N}_i=\mathbf{N}·\mathbf{1}.$$

(13)

Likewise, the bending moment follows from the balance of the angular momentum:

$$M\left(x\right)=\sum _{i=1}^n{b}_i{\int }_{z_{i-1}}^{z_i}\sigma (x,\overline{z})(\overline{z}-{\zeta }_0)\mathrm{d}\overline{z}=\sum _{i=1}^n{M}_i=\mathbf{M}·\mathbf{1}.$$

(14)

Above, the vectors $\mathbf{N}={\left\{\begin{array}{cccc}{N}_1& N_2& \cdots & N_n\end{array}\right\}}^T$ and $\mathbf{M}={\left\{\begin{array}{cccc}{M}_1& M_2& \cdots & M_n\end{array}\right\}}^T$ represent the contribution of each layer to total stress resultants.

The problem in Equations (11) and (12) can now be expressed in the setting involving the displacements and stress resultants. In the first step, the kinematics of the beam in Equation (8) is exploited to obtain:

$$\begin{array}{c}\frac{\sigma }{E_i}=-L_{\lambda ,i}^2(u_0^{\left(3\right)}-w^{\left(4\right)}(z-{\zeta }_0))+u_0^{\left(1\right)}-w^{\left(2\right)}(z-{\zeta }_0),\hfill \end{array}$$

(15)

and constraints

$$\begin{array}{c}{L}_{\lambda ,i}\left(u_0^{\left(2\right)}-w^{\left(3\right)}(z-{\zeta }_0)\right)-(u_0^{\left(1\right)}-w^{\left(2\right)}(z-{\zeta }_0))={\left.0\right|}_{x=0},\hfill \\ L_{\lambda ,i}\left(u_0^{\left(2\right)}-w^{\left(3\right)}(z-{\zeta }_0)\right)+(u_0^{\left(1\right)}-w^{\left(2\right)}(z-{\zeta }_0))={\left.0\right|}_{x=L}.\hfill \end{array}$$

(16)

The latter non-standard constraints are constitutive boundary conditions. Introducing the stress in Equation (15) into the stress resultant $N\left(x\right)$, Equation (13) yields the contribution of each layer $N_i\left(x\right)$ to the total $N\left(x\right)$:

$$\begin{array}{c}{N}_i\left(x\right)=b_i{E}_i(-L_{\lambda ,i}^2(u_0^{\left(3\right)}{h}_i-w^{\left(4\right)}{\int }_{z_{i-1}}^{z_i}(z-{\zeta }_0)\mathrm{d}\overline{z})+(u_0^{\left(1\right)}{h}_i-w^{\left(2\right)}{\int }_{z_{i-1}}^{z_i}(z-{\zeta }_0)\mathrm{d}\overline{z})).\hfill \end{array}$$

(17)

For the complete cross-section, this can be represented by means of the vectorial representation:

$$\mathbf{N}=\left\{\begin{array}{c}{N}_1\\ N_2\\ ⋮\\ N_n\end{array}\right\}=\left\{\begin{array}{c}{b}_1{E}_1\left(-L_{\lambda ,1}^2(u_0^{\left(3\right)}{h}_1-w^{\left(4\right)}{\int }_{z_0}^{z_1}(z-{\zeta }_0)\mathrm{d}\overline{z})+u_0^{\left(1\right)}{h}_1-w^{\left(2\right)}{\int }_{z_0}^{z_1}(z-{\zeta }_0)\mathrm{d}\overline{z}\right)\\ b_2{E}_2\left(-L_{\lambda ,2}^2(u_0^{\left(3\right)}{h}_2-w^{\left(4\right)}{\int }_{z_1}^{z_2}(z-{\zeta }_0)\mathrm{d}\overline{z})+u_0^{\left(1\right)}{h}_2-w^{\left(2\right)}{\int }_{z_1}^{z_2}(z-{\zeta }_0)\mathrm{d}\overline{z}\right)\\ ⋮\\ b_n{E}_n\left(-L_{\lambda ,n}^2(u_0^{\left(3\right)}{h}_n-w^{\left(4\right)}{\int }_{z_{n-1}}^{z_n}(z-{\zeta }_0)\mathrm{d}\overline{z})+u_0^{\left(1\right)}{h}_n-w^{\left(2\right)}{\int }_{z_{n-1}}^{z_n}(z-{\zeta }_0)\mathrm{d}\overline{z}\right)\end{array}\right\}.$$

(18)

Scalar multiplication by $\mathbf{1}$ gives:

$$\begin{array}{c}N\left(x\right)=-(\mathbf{A}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1}{u}_0^{\left(3\right)}+(\mathbf{S}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1}{w}^{\left(4\right)}+(\mathbf{E}\odot \mathbf{A})·\mathbf{1}{u}_0^{\left(1\right)}-(\mathbf{E}\odot \mathbf{S})·\mathbf{1}{w}^{\left(2\right)}\hfill \end{array}$$

(19)

or more compactly:

$$\begin{array}{c}N\left(x\right)=-k_{EA}^{NL}{u}_0^{\left(3\right)}+k_{ES}^{NL}{w}^{\left(4\right)}+k_{EA}{u}_0^{\left(1\right)}-k_{ES}{w}^{\left(2\right)}.\hfill \end{array}$$

(20)

A similar procedure can be applied to the constitutive boundary conditions in Equation (16) which become:

$$\begin{array}{c}({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{A}-w^{\left(3\right)}\mathbf{S}\right)-\mathbf{E}·(u_0^{\left(1\right)}\mathbf{A}-w^{\left(2\right)}\mathbf{S})={\left.0\right|}_{x=0},\hfill \\ ({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{A}-w^{\left(3\right)}\mathbf{S}\right)+\mathbf{E}·(u_0^{\left(1\right)}\mathbf{A}-w^{\left(2\right)}\mathbf{S})={\left.0\right|}_{x=L}.\hfill \end{array}$$

(21)

As for the bending moment, the same steps are followed. Starting from the normal stresses in Equation (15) that are multiplied by $(z-{\zeta }_0)$, it is obtained:

$$\begin{array}{c}M\left(x\right)=-k_{ES}^{NL}{u}_0^{\left(3\right)}+k_{EI}^{NL}{w}^{\left(4\right)}+k_{ES}{u}_0^{\left(1\right)}-k_{EI}{w}^{\left(2\right)}.\hfill \end{array}$$

(22)

The constitutive boundary conditions are:

$$\begin{array}{c}({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{S}-w^{\left(3\right)}\mathbf{I}\right)-\mathbf{E}·(u_0^{\left(1\right)}\mathbf{S}-w^{\left(2\right)}\mathbf{I})={\left.0\right|}_{x=0},\hfill \\ ({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{S}-w^{\left(3\right)}\mathbf{I}\right)+\mathbf{E}·(u_0^{\left(1\right)}\mathbf{S}-w^{\left(2\right)}\mathbf{I})={\left.0\right|}_{x=L}.\hfill \end{array}$$

(23)

This completes the nonlocal stress-driven beam relations between the displacements and stress resultants.

4. A Nonlocal Variational Setting for Beams

The governing differential equations for the beam $\mathcal{B}$ should follow from a suitable potential. Denoting the internal potential as $U_i$ and the external potential as $U_e$, the potential is then:

$$\mathsf{\Pi }\left(u_0,w\right)=U_i-U_e,$$

(24)

where

$$\begin{array}{c}{U}_i={\int }_{\mathcal{B}}{\int }_0^{\epsilon }\sigma \mathrm{d}\epsilon \mathrm{d}V,\\ U_e={\int }_L{q}_x{u}_0\mathrm{d}x+{\int }_L{q}_{z}w\mathrm{d}x+{\mathcal{N}}_0{u}_0\left(0\right)+{\mathcal{N}}_L{u}_0\left(L\right)+{\mathcal{T}}_{0}w\left(0\right)+{\mathcal{T}}_{L}w\left(L\right)-{\mathcal{M}}_0{w}^{\left(1\right)}\left(0\right)-{\mathcal{M}}_L{w}^{\left(1\right)}\left(L\right).\end{array}$$

(25)

Above, $q_z\left(x\right)$ and $q_x\left(x\right)$ represent distributed loadings in the transverse and longitudinal directions. ${\mathcal{N}}_0,{\mathcal{N}}_L$, ${\mathcal{T}}_0,{\mathcal{T}}_L$, and ${\mathcal{M}}_0,{\mathcal{M}}_L$ are the external axial forces, transversal forces, and bending moments at $x\in \{0,L\}$, respectively.

Longitudinal and transverse displacement fields follow from the principle of minimum potential energy:

$$\begin{array}{c}(u_0,w)=arg{\displaystyle \underset{u_0,w}{inf}}\mathsf{\Pi }(u_0,w),\hfill \end{array}$$

(26)

upon the enforcement of the stationarity conditions. This is pursued in the following section.

4.1. Stationarity with Respect to Axial Displacement

To obtain the axial displacements, the stationary conditions ${\delta }_{u_0}\mathsf{\Pi }=0$ is applied. With the total energy potential in Equations (24) and (13) and the strain in Equation (8), it follows that

$${\delta }_{u_0}\mathsf{\Pi }={\int }_{L}N\delta u_0^{\left(1\right)}\mathrm{d}x-{\int }_L{q}_x\delta u_0\mathrm{d}x-{\mathcal{N}}_0\delta u_0\left(0\right)-{\mathcal{N}}_L\delta u_0\left(L\right)=0.$$

(27)

Application of integration by parts gives:

$${\delta }_{u_0}\mathsf{\Pi }=N\delta u_0\left|{}_0^L\right.-{\int }_L{N}^{\left(1\right)}\delta u_0\mathrm{d}x-{\int }_L{q}_x\delta u_0\mathrm{d}x-{\mathcal{N}}_0\delta u_0\left(0\right)-{\mathcal{N}}_L\delta u_0\left(L\right)=0.$$

(28)

Assembling terms with virtual displacements $\delta u_0\left(0\right)$ and $\delta u_0\left(L\right)$ and accounting for the arbitrariness of the virtual displacements yields boundary conditions:

$$\begin{array}{cc}N\left(0\right)=-{\mathcal{N}}_0\hfill & \mathrm{or}\mathrm{prescribe}{u}_0\left(0\right),\hfill \\ N\left(L\right)={\mathcal{N}}_L\hfill & \mathrm{or}\mathrm{prescribe}{u}_0\left(L\right).\hfill \end{array}$$

(29)

Likewise, the terms with $\delta u_0$ in Equation (28) can now be assembled. Due to the arbitrariness of the virtual displacement, it is:

$$N^{\left(1\right)}+q_x=0,$$

(30)

which serves as a governing equation of the problem. The first derivative of the normal force follows from Equation (20):

$$\begin{array}{c}{N}^{\left(1\right)}=-k_{EA}^{NL}{u}_0^{\left(4\right)}+k_{ES}^{NL}{w}^{\left(5\right)}+k_{EA}{u}_0^{\left(2\right)}-k_{ES}{w}^{\left(3\right)},\hfill \end{array}$$

(31)

which now gives the final form of the equation governing the axial equilibrium from Equation (30):

$$-k_{EA}^{NL}{u}_0^{\left(4\right)}+k_{ES}^{NL}{w}^{\left(5\right)}+k_{EA}{u}_0^{\left(2\right)}-k_{ES}{w}^{\left(3\right)}+q_x=0.$$

(32)

In addition, the boundary conditions in Equation (29) can be rearranged by the application of Equation (20):

$$\begin{array}{cc}{\left.(-k_{EA}^{NL}{u}_0^{\left(3\right)}+k_{ES}^{NL}{w}^{\left(4\right)}+k_{EA}{u}_0^{\left(1\right)}-k_{ES}{w}^{\left(2\right)})\right|}_{x=0}=-{\mathcal{N}}_0\hfill & \mathrm{or}\mathrm{prescribe}{u}_0\left(0\right),\hfill \\ {\left.(-k_{EA}^{NL}{u}_0^{\left(3\right)}+k_{ES}^{NL}{w}^{\left(4\right)}+k_{EA}{u}_0^{\left(1\right)}-k_{ES}{w}^{\left(2\right)})\right|}_{x=L}={\mathcal{N}}_L\hfill & \mathrm{or}\mathrm{prescribe}{u}_0\left(L\right).\hfill \end{array}$$

(33)

These boundary conditions must be accompanied by the constitutive boundary conditions in Equation (21). In this manner, the stress resultants are eliminated from the formulation and only displacements are involved. Nevertheless, the normal force can be recovered in the post-processing phase from Equation (20) if required. For the solution procedure involving both displacements and stress resultants, see the works in [47,48,49].

4.2. Stationarity with Respect to Transverse Displacement

However, it remains to invoke the stationary condition with respect to the transverse displacement ${\delta }_w\mathsf{\Pi }=0$ to describe bending of the nonlocal beam. Following the same steps as for the longitudinal displacement, but with application of the equilibrium equation (Equation (14)), it is obtained:

$$\begin{array}{c}{\delta }_w\mathsf{\Pi }=-{\int }_{L}M\delta w^{\left(2\right)}\mathrm{d}x-{\int }_L{q}_z\delta w\mathrm{d}x-{\mathcal{T}}_0\delta w\left(0\right)-{\mathcal{T}}_L\delta w\left(L\right)+{\mathcal{M}}_0\delta w^{\left(1\right)}\left(0\right)+{\mathcal{M}}_L\delta w^{\left(1\right)}\left(L\right)=0.\hfill \end{array}$$

(34)

Integrating by parts twice provides:

$$\begin{array}{c}{\delta }_w\mathsf{\Pi }=-M\delta w^{\left(1\right)}\left|{}_0^L\right.+M^{\left(1\right)}\delta w\left|{}_0^L\right.\hfill \\ -{\int }_L{M}^{\left(2\right)}\delta w\mathrm{d}x-{\int }_L{q}_z\delta w\mathrm{d}x-{\mathcal{T}}_0\delta w\left(0\right)-{\mathcal{T}}_L\delta w\left(L\right)+{\mathcal{M}}_0\delta w^{\left(1\right)}\left(0\right)+{\mathcal{M}}_L\delta w^{\left(1\right)}\left(L\right)=0.\hfill \end{array}$$

(35)

The arbitrariness of $\delta w^{\left(1\right)}$ at $x\in \{0,L\}$ provides the first set of boundary conditions:

$$\begin{array}{cc}M\left(0\right)=-{\mathcal{M}}_0\hfill & \mathrm{or}\mathrm{prescribe}{w}^{\left(1\right)}\left(0\right),\hfill \\ M\left(L\right)={\mathcal{M}}_L\hfill & \mathrm{or}\mathrm{prescribe}{w}^{\left(1\right)}\left(L\right),\hfill \end{array}$$

(36)

while the corresponding part for $\delta w$ at $x\in \{0,L\}$ gives:

$$\begin{array}{cc}{M}^{\left(1\right)}\left(0\right)=-{\mathcal{T}}_0\hfill & \mathrm{or}\mathrm{prescribe}w\left(0\right),\hfill \\ M^{\left(1\right)}\left(L\right)={\mathcal{T}}_L\hfill & \mathrm{or}\mathrm{prescribe}w\left(L\right).\hfill \end{array}$$

(37)

The constitutive boundary conditions in Equation (23) must be respected as well.

The remaining terms in Equation (35) provide the governing differential equation for bending:

$$\begin{array}{c}-{\int }_L{M}^{\left(2\right)}\delta w\mathrm{d}x-{\int }_L{q}_z\delta w\mathrm{d}x=0,\hfill \end{array}$$

(38)

or

$$\begin{array}{c}-M^{\left(2\right)}-q_z=0.\hfill \end{array}$$

(39)

The second derivative of the bending moment is obtained from Equation (22). This provides:

$$\begin{array}{c}{M}^{\left(2\right)}=-k_{ES}^{NL}{u}_0^{\left(5\right)}+k_{EI}^{NL}{w}^{\left(6\right)}+k_{ES}{u}_0^{\left(3\right)}-k_{EI}{w}^{\left(4\right)}.\hfill \end{array}$$

(40)

With the above result at hand, the differential equation (Equation (39)) that governs the transverse displacement now becomes:

$$\begin{array}{c}-k_{ES}^{NL}{u}_0^{\left(5\right)}+k_{EI}^{NL}{w}^{\left(6\right)}+k_{ES}{u}_0^{\left(3\right)}-k_{EI}{w}^{\left(4\right)}+q_z=0,\hfill \end{array}$$

(41)

along with the constitutive boundary conditions in Equation (23). Following the same reasoning as in the longitudinal displacements, the boundary conditions in Equations (38) and (39) are rearranged by Equation (22):

$$\begin{array}{cc}{\left.(-k_{ES}^{NL}{u}_0^{\left(3\right)}+k_{EI}^{NL}{w}^{\left(4\right)}+k_{ES}{u}_0^{\left(1\right)}-k_{EI}{w}^{\left(2\right)})\right|}_{x=0}=-{\mathcal{M}}_0\hfill & \mathrm{or}\mathrm{prescribe}{w}^{\left(1\right)}\left(0\right),\hfill \\ {\left.(-k_{ES}^{NL}{u}_0^{\left(3\right)}+k_{EI}^{NL}{w}^{\left(4\right)}+k_{ES}{u}_0^{\left(1\right)}-k_{EI}{w}^{\left(2\right)})\right|}_{x=L}={\mathcal{M}}_L\hfill & \mathrm{or}\mathrm{prescribe}{w}^{\left(1\right)}\left(L\right)\hfill \end{array}$$

(42)

while the corresponding part for $\delta w$ at $x\in \{0,L\}$ gives:

$$\begin{array}{cc}{\left.(-k_{ES}^{NL}{u}_0^{\left(4\right)}+k_{EI}^{NL}{w}^{\left(5\right)}+k_{ES}{u}_0^{\left(2\right)}-k_{EI}{w}^{\left(3\right)})\right|}_{x=0}=-{\mathcal{T}}_0\hfill & \mathrm{or}\mathrm{prescribe}w\left(0\right),\hfill \\ {\left.(-k_{ES}^{NL}{u}_0^{\left(4\right)}+k_{EI}^{NL}{w}^{\left(5\right)}+k_{ES}{u}_0^{\left(2\right)}-k_{EI}{w}^{\left(3\right)})\right|}_{x=L}={\mathcal{T}}_L\hfill & \mathrm{or}\mathrm{prescribe}w\left(L\right).\hfill \end{array}$$

(43)

The bending moment can be conveniently obtained in the post-processing step from Equation (22).

4.3. Neutral Surface Position

The exact position of the neutral surface can be determined from the condition that the bending stress must disappear. To simplify further presentation, beams symmetric with respect to $x-z$ plane are considered. Such an assumption is consistent with the vast majority of practical applications. Now, suppose that the neutral surface is positioned in the layer i. Using Equation (9) as a starting point along with Equation (8) yields:

$$u_0^{\left(1\right)}-w^{\left(2\right)}(z-{\zeta }_0)={\int }_0^L{\varphi }_{\lambda ,i}(x-\xi )E_i^{-1}\sigma (\xi ,z)\mathrm{d}\xi .$$

(44)

Setting $\sigma =0$ in the latter equation and then integrating over the cross-section, with focus on the transverse displacement part, gives:

$$\begin{array}{c}\sum _{i=1}^n{b}_i{E}_i{\int }_{z_{i-1}}^{z_i}(\overline{z}-{\zeta }_0)\mathrm{d}\overline{z}=\mathbf{E}·\mathbf{S}=k_{\mathrm{ES}}=0\end{array}$$

(45)

or

$$\begin{array}{c}\sum _{i=1}^n{b}_i{E}_i{\int }_{z_{i-1}}^{z_i}\overline{z}\mathrm{d}\overline{z}-\sum _{i=1}^n{b}_i{E}_i{\zeta }_0{\int }_{z_{i-1}}^{z_i}\mathrm{d}\overline{z}=0.\end{array}$$

(46)

Simple transformation now provides the required neutral surface shift ${\zeta }_0$:

$$\begin{array}{c}{\zeta }_0=\frac{{\sum }_{i=1}^n{b}_i{E}_i(z_i^2-z_{i-1}^2)}{2(\mathbf{A}\odot \mathbf{E})·\mathbf{1}}=\frac{{\sum }_{i=1}^n{b}_i{h}_i{E}_i(z_i+z_{i-1})}{2k_{EA}}.\end{array}$$

(47)

Note that equality $k_{ES}=0$ arising from symmetry of the cross-section reduces coupling effects in the differential formulations in Equations (32) and (41) and the constitutive boundary conditions in Equations (21) and (23). The important case of complete decoupling of axial and transverse displacements can be achieved only if the term $k_{ES}^{NL}=(\mathbf{S}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1}$ vanishes (see Equation (4)). This can happen if one of the following conditions is met:

• all small-size parameters are equal, ${\lambda }_i=\lambda$, $i\in \{1,2,\dots ,n\}$; or
• symmetry about the $x-y$ plane in all material properties (Young’s modulus and small-size parameter) and geometry (width and height) of each layer exists.

In the latter case, the neutral surface shift ${\zeta }_0=0$. Depending on the particular beam configuration, a decoupled system of differential equations can be significantly easier to solve compared to the coupled problem. For readers’ convenience, frameworks for both coupled and decoupled formulation are summarized in

Table 1. Coupled nonlocal layered beams framework for beams with $x-z$ symmetry of the cross-section, $k_{ES}=0$.

Calculate stiffnesses:	$k_{EA}=(\mathbf{E}\odot \mathbf{A})·\mathbf{1}$,
	$k_{EI}=(\mathbf{E}\odot \mathbf{I})·\mathbf{1}$,
	$k_{EA}^{NL}=(\mathbf{A}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1},$
	$k_{ES}^{NL}=(\mathbf{S}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1},$
	$k_{EI}^{NL}=(\mathbf{I}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1}.$
Coupled formulation:	$-k_{EA}^{NL}{u}_0^{\left(4\right)}+k_{ES}^{NL}{w}^{\left(5\right)}+k_{EA}{u}_0^{\left(2\right)}+q_x=0$,
	$-k_{ES}^{NL}{u}_0^{\left(5\right)}+k_{EI}^{NL}{w}^{\left(6\right)}-k_{EI}{w}^{\left(4\right)}+q_z=0$.
Boundary conditions:	${\left.(-k_{EA}^{NL}{u}_0^{\left(3\right)}+k_{ES}^{NL}{w}^{\left(4\right)}+k_{EA}{u}_0^{\left(1\right)})\right|}_{x=0}=-{\mathcal{N}}_0$	$\mathrm{or}\mathrm{prescribe}{u}_0\left(0\right),$
	${\left.(-k_{EA}^{NL}{u}_0^{\left(3\right)}+k_{ES}^{NL}{w}^{\left(4\right)}+k_{EA}{u}_0^{\left(1\right)})\right|}_{x=L}={\mathcal{N}}_L$	$\mathrm{or}\mathrm{prescribe}{u}_0\left(L\right),$
	${\left.(-k_{ES}^{NL}{u}_0^{\left(3\right)}+k_{EI}^{NL}{w}^{\left(4\right)}-k_{EI}{w}^{\left(2\right)})\right|}_{x=0}=-{\mathcal{M}}_0$	$\mathrm{or}\mathrm{prescribe}{w}^{\left(1\right)}\left(0\right),$
	${\left.(-k_{ES}^{NL}{u}_0^{\left(3\right)}+k_{EI}^{NL}{w}^{\left(4\right)}-k_{EI}{w}^{\left(2\right)})\right|}_{x=L}={\mathcal{M}}_L$	$\mathrm{or}\mathrm{prescribe}{w}^{\left(1\right)}\left(L\right),$
	${\left.(-k_{ES}^{NL}{u}_0^{\left(4\right)}+k_{EI}^{NL}{w}^{\left(5\right)}-k_{EI}{w}^{\left(3\right)})\right|}_{x=0}=-{\mathcal{T}}_0$	$\mathrm{or}\mathrm{prescribe}w\left(0\right),$
	${\left.(-k_{ES}^{NL}{u}_0^{\left(4\right)}+k_{EI}^{NL}{w}^{\left(5\right)}-k_{EI}{w}^{\left(3\right)})\right|}_{x=L}={\mathcal{T}}_L$	$\mathrm{or}\mathrm{prescribe}w\left(L\right).$
Constitutive boundary	$({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{A}-w^{\left(3\right)}\mathbf{S}\right)-u_0^{\left(1\right)}\mathbf{E}·\mathbf{A}={\left.0\right|}_{x=0},$
conditions:	$({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{A}-w^{\left(3\right)}\mathbf{S}\right)+u_0^{\left(1\right)}\mathbf{E}·\mathbf{A}={\left.0\right|}_{x=L},$
	$({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{S}-w^{\left(3\right)}\mathbf{I}\right)+w^{\left(2\right)}\mathbf{E}·\mathbf{I}={\left.0\right|}_{x=0},$
	$({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{S}-w^{\left(3\right)}\mathbf{I}\right)-w^{\left(2\right)}\mathbf{E}·\mathbf{I}={\left.0\right|}_{x=L}.$

and

Table 2. Decoupled nonlocal layered beams framework for beams with material and geometric symmetry of the cross-section, $k_{ES}=0$ and $k_{ES}^{NL}=0$.

Calculate stiffnesses:	$k_{EA}=(\mathbf{E}\odot \mathbf{A})·\mathbf{1}$,
	$k_{EI}=(\mathbf{E}\odot \mathbf{I})·\mathbf{1}$,
	$k_{EA}^{NL}=(\mathbf{A}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1},$
	$k_{EI}^{NL}=(\mathbf{I}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1}.$
Axial displacements:	$-k_{EA}^{NL}{u}_0^{\left(4\right)}+k_{EA}{u}_0^{\left(2\right)}+q_x=0$
Constitutive boundary	$({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\mathbf{A}{u}_0^{\left(2\right)}-\mathbf{E}·\mathbf{A}{u}_0^{\left(1\right)}={\left.0\right|}_{x=0},$
conditions:	$({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\mathbf{A}{u}_0^{\left(2\right)}+\mathbf{E}·\mathbf{A}{u}_0^{\left(1\right)}={\left.0\right|}_{x=L},$
Boundary conditions:	${\left.(-k_{EA}^{NL}{u}_0^{\left(3\right)}+k_{EA}{u}_0^{\left(1\right)})\right|}_{x=0}=-{\mathcal{N}}_0$	$\mathrm{or}\mathrm{prescribe}{u}_0\left(0\right),$
	${\left.(-k_{EA}^{NL}{u}_0^{\left(3\right)}+k_{EA}{u}_0^{\left(1\right)})\right|}_{x=L}={\mathcal{N}}_L$	$\mathrm{or}\mathrm{prescribe}{u}_0\left(L\right).$
Transverse displacements:	$-k_{EI}^{NL}{w}^{\left(6\right)}+k_{EI}{w}^{\left(4\right)}-q_z=0$
Constitutive boundary	$-({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\mathbf{I}{w}^{\left(3\right)}+\mathbf{E}·\mathbf{I}{w}^{\left(2\right)}={\left.0\right|}_{x=0},$
conditions:	$-({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\mathbf{I}{w}^{\left(3\right)}-\mathbf{E}·\mathbf{I}{w}^{\left(2\right)}={\left.0\right|}_{x=L},$
Boundary conditions:	${\left.(k_{EI}^{NL}{w}^{\left(4\right)}-k_{EI}{w}^{\left(2\right)})\right|}_{x=0}=-{\mathcal{M}}_0$	$\mathrm{or}\mathrm{prescribe}{w}^{\left(1\right)}\left(0\right),$
	${\left.(k_{EI}^{NL}{w}^{\left(4\right)}-k_{EI}{w}^{\left(2\right)})\right|}_{x=L}={\mathcal{M}}_L$	$\mathrm{or}\mathrm{prescribe}{w}^{\left(1\right)}\left(L\right),$
	${\left.(k_{EI}^{NL}{w}^{\left(5\right)}-k_{EI}{w}^{\left(3\right)})\right|}_{x=0}=-{\mathcal{T}}_0=0$	$\mathrm{or}\mathrm{prescribe}w\left(0\right),$
	${\left.(k_{EI}^{NL}{w}^{\left(5\right)}-k_{EI}{w}^{\left(3\right)})\right|}_{x=L}={\mathcal{T}}_L=0$	$\mathrm{or}\mathrm{prescribe}w\left(L\right).$

.

5. Examples

5.1. Cantilever Beam

Consider a cantilever beam with length $L=100$ nm, loaded at the free end ($x=L$) with the longitudinal force $F_x=1$ nN and transverse force $F_z=1$ nN. The beam is composed of two layers with the following geometrical and material properties:

$$\begin{array}{c}\mathbf{b}=\left\{\begin{array}{c}1\\ 1\end{array}\right\}\mathrm{nm},\mathbf{h}=\left\{\begin{array}{c}5\\ 5\end{array}\right\}\mathrm{nm},\mathbf{\lambda }=\left\{\begin{array}{c}0.1\\ 0.2\end{array}\right\},\mathbf{E}=\left\{\begin{array}{c}{10}^6\\ 2·{10}^6\end{array}\right\}{\mathrm{nN}/\mathrm{nm}}^2.\end{array}$$

(48)

Consequently, z coordinates defining the beginning and the end of each layer are:

$$\begin{array}{c}\left\{\begin{array}{c}{z}_0\\ z_1\\ z_2\end{array}\right\}=\left\{\begin{array}{c}-5\\ 0\\ 5\end{array}\right\}\mathrm{nm}.\end{array}$$

(49)

Geometry of the beam cross-section is symmetric about $x-z$ and $x-y$ plane, but is not symmetric with respect to material properties. Therefore, $k_{ES}=0$ and

$$\begin{array}{c}{k}_{ES}^{NL}=(\mathbf{S}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1}=-\frac{b_1{b}_2{E}_1{E}_2{(h_1+h_2)}^2{L}^2({\lambda }_1^2-{\lambda }_2^2)}{4(b_1{E}_1+b_2{E}_2)}=5\times {10}^9\mathrm{nN}{\mathrm{nm}}^3.\hfill \end{array}$$

(50)

and the framework in Table 1 applies. Remaining local and nonlocal stiffnesses in Equations (3) and (4) are:

$$\begin{array}{c}{k}_{EA}=(\mathbf{E}\odot \mathbf{A})·\mathbf{1}=1.5\times {10}^7\mathrm{nN},\hfill \\ k_{EI}=(\mathbf{E}\odot \mathbf{I})·\mathbf{1}=1.15\times {10}^8\mathrm{nN}{\mathrm{nm}}^2,\hfill \end{array}$$

(51)

and

$$\begin{array}{c}{k}_{EA}^{NL}=(\mathbf{A}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1}=4.5\times {10}^9\mathrm{nN}{\mathrm{nm}}^2,\hfill \\ k_{EI}^{NL}=(\mathbf{I}\odot \mathbf{E}\odot {\mathbf{L}}_{\lambda }\odot {\mathbf{L}}_{\lambda })·\mathbf{1}=2.6\times {10}^{10}\mathrm{nN}{\mathrm{nm}}^4.\hfill \end{array}$$

(52)

The shift of the neutral surface ${\zeta }_0$ is then:

$$\begin{array}{c}{\zeta }_0=\frac{h}{4}\frac{b_2{E}_2-b_1{E}_1}{b_1{E}_1+b_2{E}_2}=0.833\mathrm{nm}.\hfill \end{array}$$

(53)

With these data, the differential formulation of the problem can be defined. According to Table 1, the governing set of differential equation to be solved is:

$$\begin{array}{c}-k_{EA}^{NL}{u}_0^{\left(4\right)}+k_{ES}^{NL}{w}^{\left(5\right)}+k_{EA}{u}_0^{\left(2\right)}=0,\hfill \\ -k_{ES}^{NL}{u}_0^{\left(5\right)}+k_{EI}^{NL}{w}^{\left(6\right)}-k_{EI}{w}^{\left(4\right)}=0.\hfill \end{array}$$

(54)

The boundary conditions are:

$$\begin{array}{c}{u}_0\left(0\right)=0,\hfill \\ {\left.(-k_{EA}^{NL}{u}_0^{\left(3\right)}+k_{ES}^{NL}{w}^{\left(4\right)}+k_{EA}{u}_0^{\left(1\right)})\right|}_{x=L}=F_x,\hfill \\ w\left(0\right)=0,\hfill \\ w^{\left(1\right)}\left(0\right)=0,\hfill \\ {\left.(-k_{ES}^{NL}{u}_0^{\left(4\right)}+k_{EI}^{NL}{w}^{\left(5\right)}-k_{EI}{w}^{\left(3\right)})\right|}_{x=0}=0,\hfill \\ {\left.(-k_{ES}^{NL}{u}_0^{\left(4\right)}+k_{EI}^{NL}{w}^{\left(5\right)}-k_{EI}{w}^{\left(3\right)})\right|}_{x=L}=F_z,\hfill \end{array}$$

(55)

while the constitutive boundary conditions are:

$$\begin{array}{c}({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{A}-w^{\left(3\right)}\mathbf{S}\right)-u_0^{\left(1\right)}\mathbf{E}·\mathbf{A}={\left.0\right|}_{x=0},\hfill \\ ({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{A}-w^{\left(3\right)}\mathbf{S}\right)+u_0^{\left(1\right)}\mathbf{E}·\mathbf{A}={\left.0\right|}_{x=L},\hfill \\ ({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{S}-w^{\left(3\right)}\mathbf{I}\right)+w^{\left(2\right)}\mathbf{E}·\mathbf{I}={\left.0\right|}_{x=0},\hfill \\ ({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\left(u_0^{\left(2\right)}\mathbf{S}-w^{\left(3\right)}\mathbf{I}\right)-w^{\left(2\right)}\mathbf{E}·\mathbf{I}={\left.0\right|}_{x=L}.\hfill \end{array}$$

(56)

The problem was solved by the aid of Wolfram Mathematica software. The specific form of the function describing the axial displacement is too long to be presented here in detail. Instead, the graphical representation is provided in Figure 2. It is well known that the local formulation exhibits linear behavior. In the nonlocal formulation marked curvature is observed providing a nonlinear distribution. Similar effects are noted elsewhere [50,51] in the presence of thermal effects. If the nonlocal parameters ${\lambda }_1,{\lambda }_2$ are allowed to vary, the obtained surface describing the maximal elongation of the beam can be graphically presented (Figure 3). Smaller values of nonlocal parameters lead toward steep surface’s gradients, but as the nonlocal parameters are increased the effect diminishes. Interestingly, for a certain choice of nonlocal parameters average axial displacements, $u_0\left(L\right)$ take negative values, although the nanobeam is loaded by the tensile axial force. This is a direct consequence of coupling terms associated to $k_{ES}^{NL}$ and the neutral surface shift. Finally, the distribution of the normal force along the beam $N=F_A=1$ nN is constant as expected (Figure 4).

The curve representing the deformed shape of the beam for given nonlocal parameters is provided in Figure 5. The expected shape is obtained. The maximal transverse displacement at the beam’s tip for different values of the nonlocal parameters ${\lambda }_1,{\lambda }_2$ is given in Figure 6. Similar conclusions about gradients as in the case of longitudinal displacement can be drawn. The bending moment is distributed linearly along the beam (Figure 4). The moment is recovered in the post-processing from Equation (22) as $M=F_z(x-L)=-100+x$ nN nm.

Statically determined problems, such as this one, show no dependence of support reactions on the nonlocal parameters (Figure 7).

5.2. Doubly Clamped Beam

The second example presents a statically indeterminate problem. A doubly clamped beam $L=500$ nm long is subjected to longitudinal constant distributed loading $q_x=1$ nN/nm and distributed transversal loading $q_z=1$ nN/nm. The beam is assembled of three layers. The first and third layers have identical geometrical and material properties. These data are provided below:

$$\begin{array}{c}\mathbf{b}=\left\{\begin{array}{c}3\\ 1\\ 3\end{array}\right\}\mathrm{nm},\mathbf{h}=\left\{\begin{array}{c}10\\ 30\\ 10\end{array}\right\}\mathrm{nm},\mathbf{\lambda }=\left\{\begin{array}{c}0.1\\ 0.2\\ 0.1\end{array}\right\},\mathbf{E}=\left\{\begin{array}{c}{10}^6\\ 2·{10}^6\\ {10}^6\end{array}\right\}{\mathrm{nN}/\mathrm{nm}}^2.\end{array}$$

(57)

Consequently, z coordinates defining the beginning and the end of each layer are:

$$\begin{array}{c}\left\{\begin{array}{c}{z}_0\\ z_1\\ z_2\\ z3\end{array}\right\}=\left\{\begin{array}{c}-25\\ -15\\ 15\\ 25\end{array}\right\}\mathrm{nm}.\end{array}$$

(58)

The beam’s cross-section is geometrically symmetric with respect to both $x-y$ and $x-z$ plane as well as with respect to material properties, thus $k_{ES}=0$ and $k_{ES}^{NL}=0$. With these data, the decoupled differential formulation of the problem can be set. Again, according to Table 2, the governing differential equation to be solved for the axial displacements is:

$$-k_{EA}^{NL}{u}_0^{\left(4\right)}+k_{EA}{u}_0^{\left(2\right)}+q_x=0.$$

(59)

The boundary conditions are:

$$\begin{array}{c}{u}_0\left(0\right)=0,\hfill \\ u_0\left(L\right)=0.\hfill \end{array}$$

(60)

The constitutive boundary conditions are:

$$\begin{array}{c}({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\mathbf{A}{u}_0^{\left(2\right)}-\mathbf{E}·\mathbf{A}{u}_0^{\left(1\right)}={\left.0\right|}_{x=0},\hfill \\ ({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\mathbf{A}{u}_0^{\left(2\right)}+\mathbf{E}·\mathbf{A}{u}_0^{\left(1\right)}={\left.0\right|}_{x=L}.\hfill \end{array}$$

(61)

As for the bending problem, the corresponding differential framework is:

$$\begin{array}{c}-k_{EI}^{NL}{w}^{\left(6\right)}+k_{EI}{w}^{\left(4\right)}+q_z=0,\hfill \end{array}$$

(62)

augmented with the following boundary conditions:

$$\begin{array}{c}w\left(0\right)=0,w^{\left(1\right)}\left(0\right)=0,\hfill \\ w\left(L\right)=0,w^{\left(1\right)}\left(L\right)=0,\hfill \end{array}$$

(63)

and constitutive boundary conditions:

$$\begin{array}{c}-({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\mathbf{I}{w}^{\left(3\right)}+\mathbf{E}·\mathbf{I}{w}^{\left(2\right)}={\left.0\right|}_{x=0},\hfill \\ -({\mathbf{L}}_{\lambda }\odot \mathbf{E})·\mathbf{I}{w}^{\left(3\right)}-\mathbf{E}·\mathbf{I}{w}^{\left(2\right)}={\left.0\right|}_{x=L}.\hfill \end{array}$$

(64)

The distribution of the axial displacement along the beam is given in Figure 8. The middle point of the beam ($x=L/2$) experiences the largest displacement, while the ends fulfill prescribed boundary conditions. Next, the movement of the middle point is analyzed for dependence on the nonlocal parameters (Figure 9). The increase in either nonlocal parameter does lead toward maximal axial displacement converging toward a limit value. The linear distribution of the normal force along the beam is defined as $N\left(x\right)=F_A-q_{x}x=250-x$ nN (Figure 10). Support reactions are obtained as $N\left(0\right)=250$ nN and $N\left(500\right)=-250$ nN.

The bell-shaped curve represents the deformed shape of the beam due to bending (Figure 11). The maximal transverse displacement takes place at the middle of the beam as expected. Influence of the nonlocal parameters ${\lambda }_1={\lambda }_3,{\lambda }_2$ is given in Figure 12. Higher gradients are obtained for lower values of the nonlocal parameter ${\lambda }_1$ compared to ${\lambda }_2$. The bending moment distribution in Equation (22) is defined by the parabolic polynomial $M\left(x\right)=-\frac{1}{2}{x}^2+250x-22,156.1$ (Figure 13), giving $M_A=M\left(0\right)=M_B=M\left(L\right)=22,156.1$ nN nm with the maximum $M=9093.85$ nN nm at $x=L/2=250$ nm. In addition, contrary to the previous example, in the statically indeterminate problems, support reactions do depend on the nonlocal parameters (Figure 14).

6. Conclusions

The following main results of the presented research are summarized below.

• The stress-driven integral approach based on Bernoulli–Euler kinematical hypotheses is extended to composite beams assembled of multiple layers, not necessarily of equal width. As demonstrated in the examples, the approach does not suffer from paradoxes present in some other formulations.
• The more standard approach that includes mixed boundary conditions, i.e., both stress resultants and prescribed displacements, is replaced by the purely kinematical framework. In this way, it is not necessary to explicitly determine support reactions in order to calculate displacements. Support reactions and stress resultant distributions are conveniently calculated in the post-processing phase.
• The example section demonstrates that in statically undetermined structural problems, reaction systems exhibit technically significant size effects which therefore have to be taken in due account in design and optimization of a wide variety of new-generation sensors and actuators.
• In the general case of layered beams, the resulting formulation exhibits coupling between axial and transverse displacements. This gives rise to unusual nonlocal phenomena, such as shortening of the nanobeam in the presence of tensile axial force. Coupling of axial and bending terms in the governing differential equations, as well as the neutral surface shift, give rise to such effects.
• Finally, as discussed in the Introduction, if the beams with larger length/thickness ratios are to be considered, one must be wary about the surface and interface effects. An extension with a specialized size-dependent model is recommended in such cases.