Stress-driven modeling of nonlocal thermoelastic behavior of nanobeams

TL;DR

This paper introduces a novel stress-driven nonlocal integral model for analyzing the thermoelastic behavior of nanobeams, overcoming limitations of traditional strain-driven models and providing exact solutions and benchmarks.

Contribution

It develops a stress-driven nonlocal thermoelastic model for nanobeams, offering an alternative to strain-driven approaches and establishing equivalence with differential equations.

Findings

Model successfully avoids strain-driven issues

Exact solutions and new benchmarks provided

Demonstrates improved accuracy in nano- and microbeam analysis

Abstract

A consistent stress-driven nonlocal integral model for nonisothermal structural analysis of elastic nano- and microbeams is proposed. Most nonlocal models of literature are strain-driven and it was shown that such approaches can lead toward a number of difficulties. Following recent contributions within the isothermal setting, the developed model abandons the classical strain-driven methodology in favour of the modern stress-driven elasticity theory by G. Romano and R. Barretta. This effectively circumvents issues associated with strain-driven formulations. The new thermoelastic nonlocal integral model is proven to be equivalent to an adequate set of differential equations, accompanied by higher-order constitutive boundary conditions, when the special Helmholtz averaging kernel is adopted in the convolution. The example section provides several applications, thus enabling insight into performance of the formulation. Exact nonlocal solutions are established, detecting also new benchmarks for thermoelastic numerical analyses.

Equations

$$u_x(x,z)=u(x,z)=u_0(x)+\phi (x)z,u_z(x)=w(x),u_y=0,$$

$$u_0(x)=\frac{1}{A}{\int }_{\Omega }u(x,z)\mathrm{d}A.$$

$$0={\gamma }_{\mathrm{xz}}(x)={\partial }_{x}w(x)+{\partial }_{z}u(x,z)=w^{(1)}(x)+\phi (x),$$

$$\epsilon (x,z)={\partial }_{x}u(x,z)=u_0^{(1)}(x)+{\phi }^{(1)}(x)z=u_0^{(1)}(x)-w^{(2)}(x)z.$$

$${\epsilon }_0=u_0^{(1)}.$$

$$w^{(2)}=\chi .$$

$$\epsilon (x,z)={\epsilon }_0(x)-z\chi (x).$$

$$\epsilon ={\epsilon }_{\mathrm{th}}+{\epsilon }_{\mathrm{el}}.$$

$$\alpha \Delta \theta (x,z)+\frac{\sigma (x,z)}{E}={\epsilon }_0(x)-z\chi (x),$$

$$\begin{array}[]{l}\displaystyle\int_{\Omega}\sigma dA=0\quad\rightarrow\quad\varepsilon_{0}=\frac{1}{AE}\displaystyle\int_{\Omega}\alpha E\Delta\theta dA,\vspace{6pt}\cr\displaystyle\int_{\Omega}\sigma zdA=0\quad\rightarrow\quad\chi=-\frac{1}{I_{\mathrm{y}}E}\displaystyle\int_{\Omega}\alpha E\Delta\theta zdA.\end{array}$$

$$\sigma =-\alpha \Delta \theta E+\frac{{\int }_{\Omega }\alpha \Delta \theta E\mathrm{d}A}{A}+\frac{{\int }_{\Omega }\alpha E\Delta \theta z\mathrm{d}A}{I_{\mathrm{y}}}z,$$

$$\sigma =-\alpha \Delta \theta E+\frac{N+{\int }_{\Omega }\alpha \Delta \theta E\mathrm{d}A}{A}+\frac{M+{\int }_{\Omega }\alpha E\Delta \theta z\mathrm{d}A}{I_{\mathrm{y}}}z,$$

$$w^{(2)}=-\frac{M+{\int }_{\Omega }\alpha E\Delta \theta z\mathrm{d}A}{E{I}_{\mathrm{y}}},$$

$${\epsilon }_{\mathrm{el}}(x,z)={\int }_0^L\varphi (x-\xi )E^{-1}\sigma (\xi ,z)\mathrm{d}\xi ,$$

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

$$\alpha \Delta \theta (x,z)+{\int }_0^L\varphi (x-\xi )E^{-1}\sigma (\xi ,z)\mathrm{d}\xi ={\epsilon }_0(x)-z\chi (x),$$

$${\int }_{\Omega }\alpha \Delta \theta (x,z)\mathrm{d}A+{\int }_{\Omega }{\int }_0^L\varphi (x-\xi )E^{-1}\sigma (\xi ,z)\mathrm{d}\xi \mathrm{d}A={\int }_{\Omega }{\epsilon }_0(x)\mathrm{d}A-{\int }_{\Omega }z\chi (x)\mathrm{d}A.$$

$${\epsilon }_0(x)=\frac{1}{A}\left\{{\int }_{\Omega }\alpha \Delta \theta (x,z)\mathrm{d}A+{\int }_0^L\varphi (x-\xi )E^{-1}N(\xi )\mathrm{d}\xi \right\}.$$

$$\chi (x)=-\frac{1}{I_{\mathrm{y}}}\left\{{\int }_{\Omega }\alpha \Delta \theta (x,z)z\mathrm{d}A+{\int }_0^L\varphi (x-\xi )E^{-1}M(\xi )\mathrm{d}\xi \right\}.$$

$${\Phi }_{\mathrm{N}}(x)=\frac{1}{A}{\int }_{\Omega }\alpha \Delta \theta (x,z)\mathrm{d}A,{\Phi }_{\mathrm{M}}(x)=\frac{1}{I_{\mathrm{y}}}{\int }_{\Omega }\alpha \Delta \theta (x,z)z\mathrm{d}A.$$

$$\begin{array}[]{l}\displaystyle u_{0}^{(1)}-\varPhi_{\mathrm{N}}=\int_{0}^{L}\phi({x-\xi)}\frac{N(\xi)}{EA}\mathrm{d}\xi,\vspace{6pt}\cr\displaystyle-w^{(2)}-\varPhi_{\mathrm{M}}=\int_{0}^{L}\phi({x-\xi)}\frac{M(\xi)}{EI_{\mathrm{y}}}\mathrm{d}\xi.\end{array}$$

$$\frac{N}{L_{\lambda }^{2}EA}={\Phi }_{\mathrm{N}}^{(2)}-u_0^{(3)}+\frac{1}{L_{\lambda }^2}\left(u_0^{(1)}-{\Phi }_{\mathrm{N}}\right).$$

$$\begin{array}[]{l}\displaystyle u_{0}^{(2)}(0)-\varPhi_{\mathrm{N}}^{(1)}(0)-\frac{1}{L_{\lambda}}\left(u_{0}^{(1)}(0)-\varPhi_{\mathrm{N}}(0)\right)=0,\vspace{6pt}\cr\displaystyle u_{0}^{(2)}(L)-\varPhi_{\mathrm{N}}^{(1)}(L)+\frac{1}{L_{\lambda}}\left(u_{0}^{(1)}(L)-\varPhi_{\mathrm{N}}(L)\right)=0.\end{array}$$

$$\frac{M}{L_{\lambda }^{2}E{I}_{\mathrm{y}}}=w^{(4)}+{\Phi }_{\mathrm{M}}^{(2)}+\frac{1}{L_{\lambda }^2}(-w^{(2)}-{\Phi }_{\mathrm{M}}),$$

$$\begin{array}[]{l}\displaystyle-w^{(3)}(0)-\varPhi_{\mathrm{M}}^{(1)}(0)-\frac{1}{L_{\lambda}}\left(-w^{(2)}(0)-\varPhi_{\mathrm{M}}(0)\right)=0,\vspace{6pt}\cr\displaystyle-w^{(3)}(L)-\varPhi_{\mathrm{M}}^{(1)}(L)+\frac{1}{L_{\lambda}}\left(-w^{(2)}(L)-\varPhi_{\mathrm{M}}(L)\right)=0.\\ \end{array}$$

$$f(x)={\int }_{x_1}^{x_2}\varphi (x-y)g(y)\mathrm{d}y,\mathrm{where}\varphi (x-y)=\frac{1}{2L_{\lambda }}\exp (-\frac{\left|x-y\right|}{L_{\lambda }}),$$

$$\begin{array}[]{l}f_{1}(x)=\displaystyle\int_{x_{1}}^{x}\phi(x-y)g(y)\mathrm{d}y=\displaystyle\int_{x_{1}}^{x}\frac{1}{2L_{\lambda}}\exp(\frac{y-x}{L_{\lambda}})g(y)\mathrm{d}y,\vspace{6pt}\cr f_{2}(x)=\displaystyle\int_{x}^{x_{2}}\phi(x-y)g(y)\mathrm{d}y=\displaystyle\int_{x}^{x_{2}}\frac{1}{2L_{\lambda}}\exp(\frac{x-y}{L_{\lambda}})g(y)\mathrm{d}y,\\ \end{array}$$

$$\begin{array}[]{l}f_{2}(x_{1})=f(x_{1}),\quad f_{1}(x_{1})=0,\\ f_{1}(x_{2})=f(x_{2}),\quad f_{2}(x_{2})=0.\end{array}$$

$$\begin{array}[]{l}\displaystyle f_{1}^{\prime}(x)=\frac{1}{2L_{\lambda}}g(x)-\frac{1}{L_{\lambda}}\displaystyle\int_{x_{1}}^{x}\frac{1}{2L_{\lambda}}\exp(\frac{y-x}{L_{\lambda}})g(y)\mathrm{d}y=\frac{1}{2L_{\lambda}}g(x)-\frac{1}{L_{\lambda}}f_{1}(x),\vspace{6pt}\cr\displaystyle f_{2}^{\prime}(x)=-\frac{1}{2L_{\lambda}}g(x)+\frac{1}{L_{\lambda}}\displaystyle\int_{x}^{x_{2}}\frac{1}{2L_{\lambda}}\exp(\frac{x-y}{L_{\lambda}})g(y)\mathrm{d}y=-\frac{1}{2L_{\lambda}}g(x)+\frac{1}{L_{\lambda}}f_{2}(x).\\ \end{array}$$

$$f^{\prime }(x)=f_1^{\prime }(x)+f_2^{\prime }(x)=\frac{1}{L_{\lambda }}(f_2(x)-f_1(x)).$$

Full text

††journal: International Jounral of Engineering Sciences

1 Introduction

Research breakthroughs in nanotechnology over recent years have caused an increased interest in the mechanical behaviour of structures at nanoscale as well. Most commonly, sensors that measure forces or displacements are analysed. However, the mechanics of such structures fails outside usual macroscale mechanical principles. Origins of such behaviour are manifold. The structures at the nanoscale have discrete nature manifested in the form of atoms and interactions between atoms so that continuum mechanics can have limited success in describing discrete nanostructures. Furthermore, forces that are completely irrelevant at the macroscale, van der Waals forces for example, can dominate the nanoscale behaviour. As a consequence, size effects start to appear. The research community has been very active lately trying to capture this behaviour by accounting for the nonlocal nature of the phenomenon. Although such a claim can be made for the problems involving isothermal deformation processes, when the nonisothermal problems are concerned the results are not so numerous. For a short review of the existing nonisothermal models, see [8] for statical and [36] for dynamical problems.

The existing methods are almost exclusively based on gradient methods and are applied to a variety of different problems, see [5, 30, 32, 21, 29, 31] for a start. As a cornerstone of most approaches, the assumption that the strain gradient also contributes to the stress state in a point is utilized [18]. To tune-up theoretical models to experimentally observed behaviour, the existence of the so-called nonlocal (small-size) constitutive parameter is postulated. However, the literature survey in [4] shows surprisingly low number of results on determination of the nonlocal parameter obtained either by means of an experimental or a theoretical procedure. The latter paper also shows that nonlocal behaviour of the nanostructures do exists even in the case of harmonic interatomic potential and may be attributed to the discrete nature of the structure. Moreover, it demonstrates that the gradient methods have difficulties matching simulated bending of carbon nanotubes.

A gradient method, widely adopted to describe size-dependent phenomena in nanostructures, is based on Eringen differential model (EDM) associated with the strain-driven nonlocal integral theory conceived in [12]. As shown more than a decade ago [22], although nanosensors are usually designed as a cantilever nanobeam with a tip force and nonlocal effects are readily experimentally observed, EDM is not adequate to assess size effects. Elastic responses associated with EDM are technically unacceptable, as discussed in [10, 9, 14, 11, 17]. Strain-driven nonlocal integral theory, introduced and successfully adopted by Eringen to study (in unbounded domains) screw dislocations and surface waves, is inapplicable to Structural Mechanics [24].

This conclusion is due to the fact that the elastostatic problem of a continuous structure defined on a bounded domain, formulated by the strain-driven nonlocal integral model, admits no solution for all static schemes of engineering interest. Nonlocal (strain-driven) stress fields in bounded structural domains are indeed not included in the affine manifold of equilibrated stresses fields, as recognized in literature [34, 2, 13, 33, 38, 15, 19, 39, 37, 35, 16] on the basis of the original contribution in [28]. All difficulties can be overcome by resorting to the innovative stress-driven nonlocal integral theory recently proposed by G. Romano and R. Barretta [25]. According to the stress-driven approach, the nonlocal elastic strain field is the convolution between the stress field and a suitable averaging kernel. Properties and merits of the stress-driven strategy in comparison with strain-driven formulations can be found in [26, 27]. Transverse free vibrations of Bernoulli-Euler nanobeams are investigated in [1] by stress-driven integral approach.

The research at hand aims to provide a well-posed nonlocal integral formulation of nanobeam mechanics in nonisothermal regime. It is assumed that the nanobeam can be described by Bernoulli-Euler kinematics. This kind of model has not been previously addressed in literature. The method presents an extension of the contributions presented in [25] and therefore can be categorized as a stress-driven model.

2 Kinematics of Nonisothermal Bernoulli-Euler Beams

This section will introduce the notation and provide well-known governing equations that will serve as a starting point for the nonlocal integral formulation introduced in §4. In the subsequent analysis, a straight nanobeam made of the material with the coefficient of thermal expansion $\alpha$ is considered. The nanobeam’s cross-section $\Omega$ is assumed to lay in the $y-z$ plane while the longitudinal axis is denoted with $x$, Fig. 1. The longitudinal axis $x$ is assumed to pass through centroids of the cross-section. Due to the Bernoulli-Euler assumption, only normal stresses directed along $x$ axis exist and these will be denoted by $\sigma$. To simplify elaborations and not to lose focus on the central problem, bending in $x-z$ plane will be considered. Finally, the axes $y$ and $z$ are assumed to be the principal axes of inertia of the cross-section, thus giving $I_{\mathrm{yz}}=0$. As usual, $A$ is the area and $I_{\mathrm{y}}$ is the second moment of area of the cross-section, respectively.

The displacement $\mathbf{u}=\left\{u_{x},u_{y},u_{z}\right\}$ of any point of the beam can be evaluated as:

[TABLE]

where $\varphi(x)$ denotes the angle of rotation of the cross-section. As the above equation indicates, the longitudinal displacement $u$ is composed of two parts. The first part $u_{0}$ represents the average displacement of the cross section, defined as:

[TABLE]

When introduced into the Eq.(2), the second part of displacement field $\varphi(x)z$ in Eq. (1) gives $\varphi(x)\int_{\Omega}z\mathrm{d}A=\varphi(x)\;S_{z}$. If the origin of the axis $z$ is positioned on the symmetry line passing through the centroid of the cross-section, the first moment of area $S_{z}$ vanishes providing the result in Eq. (2).

Due to the Bernoulli-Euler assumption, the cross-section remains plane. This effectively enforces vanishing shear strains, thus providing a link between the derivatives of the longitudinal and the transversal displacement as:

[TABLE]

where Eq.(1)1 was utilized. The apex (n) denotes the $n$-th derivative with respect to the longitudinal coordinate $x$. With the link between the transversal displacement and the rotation established, the strains in a point can be now conveniently obtained by differentiating the axial displacement field with respect to the longitudinal coordinate as:

[TABLE]

The above introduced kinematic framework can account for both isothermal and nonisothermal class of problems.

Note that the normal strain can be additively separated into an axial and a bending part. The axial part of the normal strain can be naturally introduced as the centroidal normal strain by setting $z=0$ in Eq. (4), providing:

[TABLE]

Likewise, a standard procedure shows that the curvature $\chi$ of the deformed axis of the beam is related to the transverse displacement $w$ as:

[TABLE]

This allows definition of the bending part of the normal strain as $z\;\chi(x)$. Hence, the normal strain can be represented as:

[TABLE]

3 Stresses and Deflections in Local Nonisothermal Bernoulli-Euler Beams

As this point, calculation of stresses in the local Bernoulli-Euler nonisothermal beam should be summarized. The case in which the beam is free from external mechanical loads, but subjected to the non-homogeneous temperature field variation $\Delta\theta(x,z)$ is considered. In such type of problems stresses can arise as well. The first step involves separation of the normal strain into a part that describes thermal dilatation due to free expansion (or contraction) of a beam $\varepsilon_{\mathrm{th}}(x,z)=\alpha\Delta\theta(x,z)$ and the one caused by the elastic effects $\varepsilon_{\mathrm{el}}(x,z)=\sigma(x,z)/E$:

[TABLE]

It should be emphasized that only elastic (mechanical) part of the strain is used to calculate stresses; free thermal elongation does not contribute to stress evolution. Having in mind kinematics of the Bernoulli-Euler beam as described in §2, the sum of thermal and elastic strain should equal axial and bending deformation, Eq. (7):

[TABLE]

so the normal stress is obtained as $\sigma(x,z)=-\alpha E\Delta\theta(x,z)+\varepsilon_{0}(x)E-zE\;\chi(x)$. In the absence of external mechanical loading and upon introduction of the normal stress, equilibrium equations provide the axial strain and the curvature:

[TABLE]

With known $\varepsilon_{0}$ and $\chi$, the normal stress is rewritten as:

[TABLE]

see classical textbooks on the subject [6, 20] for the in-depth discussion. On the other side, if the beam is subjected to simultaneous action of mechanical and thermal loads, the stresses are obtained by the simple superposition procedure:

[TABLE]

where $N$ is the the axial force and $M$ is the bending moment in the corresponding cross-section. The two integrals $\int_{\Omega}\alpha\Delta\theta E\mathrm{d}A$ and $\int_{\Omega}\alpha E\Delta\theta z\mathrm{d}A$ are usually named as the thermal axial force $N_{\mathrm{th}}$ and the thermal bending moment $M_{\mathrm{th}}$.

The deflection line can be obtained using Eq. (6) and results Eq. (10)2, thus providing the differential equation:

[TABLE]

accompanied by a suitable set of boundary conditions.

4 Stress-driven nonlocal thermoelasticity

The framework described so far should be now extended to account for nonlocality present at the nanoscale level. The first step again involves separation of the strain into a part that describes dilatation due to free expansion (or contraction) of a beam $\varepsilon_{\mathrm{th}}(x,z)=\alpha\Delta\theta(x,z)$ and the one caused by the stresses $\varepsilon_{\mathrm{el}}$, Eq. (8). In contrast to the earlier works on the nonlocal beam mechanics cited in the Introduction, the present approach is stress-driven. Accordingly, the elastic strains originating from stresses are assumed to be of the nonlocal nature and defined as a convolution [26]:

[TABLE]

where the kernel function $\phi({x)}$ is

[TABLE]

The characteristic length is defined as $L_{\lambda}=\lambda L>0$, i.e. as the product of the nonlocal parameter $\lambda>0$ and the beam’s length $L$. Due to the exponential nature of the kernel function, the nonlocal effects in the vicinity of the corresponding point have significantly more influence on the strain than points situated at larger distance. Note that above assumption accounts only for nonlocality in the longitudinal direction. Variation of stresses in the transverse directions does not contribute to the nonlocal normal strain in the present formulation.

Since the decomposition Eq. (7) remains valid, the normal strain provides the equality:

[TABLE]

where the thermal strain and the nonlocal elastic strain definition were used. If the latter equation is integrated over the cross section, it follows:

[TABLE]

The proposition that the cross-section axes are central and principal axes imply that the second term on the right hand side vanishes. In the second term on the left hand side, the equilibrium equation $\int_{\Omega}\sigma\mathrm{d}A=N$ can be employed. This yields:

[TABLE]

In the same manner, multiplication of Eq. (16) by the transverse coordinate $z$, the equilibrium equation $\int_{\Omega}\sigma z\mathrm{d}A=M$ and subsequent integration over the cross section will provide the curvature as:

[TABLE]

Obviously, the equilibrium conditions are explicitly enforced in the formulation. In order to simplify notation, at this point we introduce the following functions:

[TABLE]

Now, accounting for Eq. (5, 6), a decoupled system of ordinary differential equations governing nonlocal mechanical behaviour of the nanobeam is obtained:

[TABLE]

The solutions of the above integrals are provided in [23]. For the average axial displacement this is:

[TABLE]

This differential equation must be augmented with the constitutive boundary conditions:

[TABLE]

For the proof of the above integral to differential transformation, see Remark 1 below for a detailed explanation. Even in the isothermal problems, the boundary conditions are of the constitutive character since they include a material property - the nonlocal parameter $L_{\lambda}$. In the nonisothermal problems, an additional material property is included in $\varPhi_{\mathrm{N}}$ and $\varPhi_{\mathrm{M}}$, the coefficient of heat expansion $\alpha$.

The transverse displacements can be dealt with in the same manner:

[TABLE]

with the boundary conditions:

[TABLE]

If the problem is statically indeterminate, the set of differential equations Eqs.(22, 24) and boundary conditions Eqs.(23, 25) can be solved upon providing additional kinematic or static boundary conditions used in the standard beam theory.

Remark 1. The equivalence of Eq. (21)1 and Eqs. (22, 23) is obtained by utilizing results [23]. However, the procedure in [23] demonstrates only necessity of existence of these conditions. To provide uniqueness of the solution, the proof in [25] in now generalized and applied to the problem at hand. In particular, consider an integral function:

[TABLE]

where $L_{\lambda}$ is a constant. In the first step, the above integral is additively separated into two parts:

[TABLE]

i.e., $f(x)=f_{1}(x)+f_{2}(x)$. In that way, the modulus in the kernel function $\phi(x-y)$ is removed. This implies following values of these functions at the boundary:

[TABLE]

By virtue of the Leibniz integral rule, differentiation of Eq. (27) with respect to the longitudinal coordinate provides:

[TABLE]

Thus, the first derivative $f^{\prime}(x)$ is then:

[TABLE]

In the same manner, second derivatives are:

[TABLE]

Now, $f^{\prime\prime}(x)=f_{1}^{\prime\prime}(x)+f_{2}^{\prime\prime}(x)$ gives:

[TABLE]

The necessary boundary conditions are obtained from values of functions at the boundary, Eqs. (28) and the first derivative Eq. (30) as:

[TABLE]

To prove the uniqueness of the solution, the solution of the homogeneous part of Eq. (32) is sought:

[TABLE]

The general form of the solution to the above differential equation is:

[TABLE]

Enforcement of the boundary conditions Eq. (33) results in:

[TABLE]

Determination of constants $C_{1}$ and $C_{2}$ from the above system gives trivial solution. Now, if there are two solutions to the non-homogeneous problem Eq.  (32) and they are denoted as $F_{1}(X)$ and $F_{2}(x)$, then the difference of two solutions must give solution of the homogeneous problem:

[TABLE]

what upon enforcement of trivial solution ensures:

[TABLE]

i.e. that the initial assumption about two different solutions was false and that only one solution exists. This proves uniqueness of the solution.

Remark 2. After determination of the displacement fields, calculation of stresses is straightforward. Namely, Eq. (16) can be transformed by introducing kinematic constraints Eq. (5, 6) into:

[TABLE]

The stress field now follow as

[TABLE]

while the accompanying constitutive boundary conditions of type Eq. (23, 25) were already accounted for in the evaluation of displacement fields.

In the absence of external forces, the thermal stresses can be alternatively evaluated from the convolution Eq. (16), along with Eqs. (18, 19) and definitions of Eq. (20) by setting $N=0,M=0$:

[TABLE]

Thermal stresses are now:

[TABLE]

It is easy to verify that if the product $\alpha\Delta\theta$ is constant, thermal stresses do not arise. The similar result follows in some other situations. For instance, when the temperature field is independent of the transverse coordinate $z$ and simultaneously linearly dependent on the longitudinal coordinate $x$ (vice-versa is also true).

5 Examples

To illustrate behaviour of the proposed methodology, several examples are considered. Introductory examples are simpler providing an interpretation of familiar topics from the standard courses of mechanics. Subsequently, more complex examples are given. The examples are solved by the aid of Wolfram Mathematica software.

5.1 Simultaneous action of uniform temperature and axial force

A homogeneous beam of a symmetric cross-section of area $A$ and Young’s modulus $E$ is subjected to a uniform temperature change $\Delta\theta$. At $x=0$ beam’s expansion is constrained, $u_{0}=0$. Otherwise the beam can expand freely. It is additionally loaded with the tensile longitudinal force $P$ at the beam’s other end, $x=L$. It is necessary to determine the extension and stresses in the beam.

Due to symmetricity of the cross section, first moment of the area $S_{y}$ vanishes. Thus, functions $\varPhi_{\mathrm{N}}(x)$ and $\varPhi_{\mathrm{M}}(x)$, Eqs.(20) are:

[TABLE]

Since $\varPhi_{\mathrm{M}}(x)=0$, the transverse displacement $w$ equals zero for the whole domain. On the other side, the elastic (mechanical) axial force field is non-vanishing, homogeneous and equal to $N=P$. Equilibrium ensures that the internal force is not a function of axial coordinate, $N\neq f(x)$, so the axial displacement Eq.(21) is:

[TABLE]

where Eq. (43)1 was used. The solution follows from Eq.(22) :

[TABLE]

where $(\alpha\Delta\theta)^{(2)}=0$ due to constant heat expansion coefficient and temperature. This can be augmented with the constitutive boundary conditions Eq. (23):

[TABLE]

Integration of the above problem yields:

[TABLE]

The first and the second part represent the well-known local result while the third part accounts for the nonlocality. To provide a graphical insight into distribution of the axial displacement, values $P=1,L=1,EA=1,\alpha=0.1,\Delta\theta=0.1$ are chosen. Results for different values of $\lambda$ are given in Fig. 2. A slight deviation from the linear distribution is noticeable at the beam’s ends. The continuous reduction of the axial displacement under constant loading with the increase of the nonlocal parameter implies that increase of the nonlocal parameter effectively increases the beam stiffness.

The nonlocal strain is obtained by differentiation as:

[TABLE]

Straightforward application of the Remark 2 enables simple calculation of stress from the equilibrium equation as $\sigma=P/A$. So, although the obtained stress is equal to the local solution, the resulting strain and elongation are distinctly nonlocal.

Finally, the resulting axial displacement field could be solved by resorting to analysis of two decoupled problems. In particular, one can solve for thermal displacements in the first step, while the second step involves solution due to action of mechanical loading. The total solution is obtained as a linear combination of these results.

5.2 Uniformly heated doubly clamped beam

In this example, a doubly clamped beam of similar geometrical and material properties like in the previous example is considered. The beam is clamped at both ends so it cannot expand. The total axial deformation is

[TABLE]

where the second - elastic part reflects the nonlocality caused by the reactions $P$ in supports. Also, the axial force does not change along the beam, so $N=P=\mathrm{const}$. Introducing nonlocality for the elastic part in the same manner as in Example 5.1 and exploiting Eq. (49), gives:

[TABLE]

Considering Eq.(22), the problem can be now stated as:

[TABLE]

augmented with two boundary conditions:

[TABLE]

Additionally, since the beam is clamped it must be $u_{0}(0)=0$. The differential equation is of the third order, so only three boundary conditions can be accounted for at this point. The remaining one, $u_{0}(L)=0$, will be introduced at the later stage. Because this problem has the identical structure as the one in Example 5.1, resulting displacement distribution is identical, cf. Eq.(47):

[TABLE]

with the difference that this equation contains unknown constant - reaction $P$. The reaction can be determined by exploiting condition $u_{0}(L)=0$ in Eq. (53), giving:

[TABLE]

Introducing the reaction $P$ into Eq. (53) gives the final distribution of the axial displacement:

[TABLE]

Contrary to the statically determine Example 5.1, this solution involves coupled thermal and elastic part that cannot be separated into two decoupled problems in the manner suggested in the former example.

To present results graphically, the same parameters as in Example 5.1 were used. Obtained total axial displacements $u_{0}$ are given in Fig. 3 for a set of nonlocal parameters. Since the local formulation will lead toward homogeneous solution $u_{0}(x)=0$ (simulated here with $\lambda=0.0001$), the effect of nonlocality is evident through the non-homogeneous displacement field. Inspection of Fig. 3 reveals that the nanobeam appears to initially becomes more flexible manifested by larger axial displacements, but after some threshold value of the nonlocal parameter, it becomes more rigid. It is also interesting to note that the nonlocal solution for the increasing value of the nonlocal parameter ($\lambda=10$) converge to the local solution ($\lambda=0$). The three-dimensional representation, Fig. 4, shows that the lowest turning point is slightly below $\lambda=0.2$. The starting softening behaviour is due to the fact that the local axial displacements are vanishing, so that for small values of $\lambda$ the nonlocal axial displacements are different from zero but as $\lambda$ tends to infinity the stress-driven total deformation goes to zero. A starting soft behaviour must be associated with any nonlocal model (strain-driven or stress-driven formulations) since local displacement (in this example) are vanishing. Similar behaviour was previously noticed in [26] in the bending of a beam. An important conclusion observed from the present example is that the terms "stiffer" and "softer" (widely used in literature) with respect to the nonlocal parameter are often misleading and often meaningless. This observation is motivated from the fact that the nonlocal total deformation is a field, so that the sign changes are often unpredictable.

Results presented in the previous example showed that increase in the nonlocal parameter will increase the nanobeam stiffness. However, in this example, the nonlocal elastic shortening at the nanobeam ends must match the thermal elongation that has the local character. This means the support reactions cannot remain constant but have to increase as well in order to compensate for increased stiffness due to $\lambda$. Fig. 5 clearly illustrates this dependency on the nonlocal parameter. Therefore, due to presence of additional kinematical constraints that relate nonlocal kinematical quantity (mechanical shortening due to support reactions) to a local kinematical quantity (thermal elongation) in statically indeterminate problems, one inevitably obtains reactions that depend on the constitutive behaviour. The same constitutive dependence can be observed even in the local case, where support reactions depend on the coefficient of heat expansion.

Note that curve $P$ vs. $\lambda$ is initially slightly curved and later almost linear, Fig. 5. This can be explained by a careful analysis of Eq. (53). In particular, the first and the second term are linear, while the third is nonlinear. Hence, the third term is responsible for reduction/increase of stiffness. For the smaller values of $\lambda$, the third term has some influence, visible as the initial nonlinearity in Fig. 5. As $\lambda$ increases, exponential terms change only slightly so the third term becomes almost a constant. As a consequence, the sum of the second and the third term tends to be a linear form for larger $\lambda$. A comparison of behaviour of the third term for a constant value of $P$ (i.e. Example 5.1) and for a variable one from the present example is provided in Fig. 6. Such behaviour was not noticed in [3] in the isothermal statically indeterimine problems, where monotonic increase in the nanobeam stiffness with the increase in the nonlocal parameter was observed. The principal difference between these two cases is the fact that in the isothermal statically indetermine problems only prescribed displacement in a point (usually equal to zero) is given. In the present case, a part of solution (thermal elongation) is known in advance for the whole domain, so this could be considered as an additional constraint on the allowable forms of axial displacement.

The normal strain $\varepsilon(x)$ distribution is given in Fig. 7. Since the local formulation results in the homogeneous normal strain field, in the present case nonlocal behaviour for $\lambda=0.1$ is clearly visible once more. Inset in the top right corner presents the homogeneous stress field as obtained from this formulation.

5.3 Bending of a nonlocal beam due to nonuniform temperature field in the transverse direction

In this example, the identical problem as in [7] in considered, thus providing a basis for comparison of the gradient and the integral approach. The cantilever beam’s length is $L$, while the rectangular cross-section is defined by the height $h$ and the width $b$. The temperature field is given as $\theta(z)=\theta_{0}+a_{1}e^{a_{2}(z+\frac{h}{2})}$, where $a_{1}$ and $a_{2}$ are constants defining distribution of temperature. The beam is not loaded by a mechanical load in any way. Thus, due to its nonlinear temperature dependency on the transverse coordinate, the nanobeam will exhibit in-plane bending. It will be assumed that the initial temperature was homogeneous and equal to zero, $\theta_{0}=0$.

Functions $\varPhi_{\mathrm{N}}$ and $\varPhi_{\mathrm{M}}$, Eq. (20), are obtained from the provided temperature field as:

[TABLE]

Since the temperature field does not depend on the longitudinal coordinate, both functions are reduced to constants and their derivatives will vanish. Also, due to absence of external mechanical loads $M=0$. The transverse displacement is the solution of differential equation Eq. (24):

[TABLE]

augmented with the boundary conditions, Eq. (25):

[TABLE]

The solution of the above problem is:

[TABLE]

what is identical to the solution based on the gradient formulation, cf. [7]. The same conclusion can be stated for the rotation:

[TABLE]

The average axial displacement follows from the differential equation, Eq. (22):

[TABLE]

with the boundary conditions:

[TABLE]

The solution is:

[TABLE]

what again corresponds to the gradient solution. Obviously, the integral formulation does not show nonlocal behaviour as well, due to particular temperature distribution and absence of external mechanical loading.

Stress can be conveniently determined by application of Eq. (42). Since the temperature field is independent of the longitudinal coordinate, the expression is reduced to:

[TABLE]

This finally yields:

[TABLE]

As expected, the stresses are not affected by nonlocality and should be obtained from the local form, Eq. (12) as well. The thermal axial force and the thermal bending moment are evaluated as:

[TABLE]

The above results introduced in Eq. (12) provides desired verification. Note that stress dependency on the nonlocal parameter could be obtained if $\partial_{x}^{2}\varPhi_{\mathrm{N}}$ and $\partial_{x}^{2}\varPhi_{\mathrm{M}}$ does not vanish, Eq. (42).

5.4 Nonuniform heating of a doubly clamped beam

In this example a doubly clamped beam $L=100$ nm long is considered. The cross-section is situated in $y-z$ plane and is assumed to be of a square shape with height $h=1$ nm and width $b=1$ nm. The Young’s modulus is $E=1$ nN/nm2. The coefficient of thermal expansion is $\alpha=1$ 1/K. The temperature increase varies along the beam as:

[TABLE]

where $a_{0}=1$ K/nm2. Consequently, temperature is a function of both longitudinal coordinate $x$ and the transverse coordinate $z$. Note that the plane $z=-h/2$ does not expand due to temperature variation. As a consequence, $\varPhi_{\mathrm{N}}$ and $\varPhi_{\mathrm{M}}$ are now functions of $x$:

[TABLE]

and their first and second derivative do not vanish. At the coordinate system’s origin $x=0$, both functions and their derivatives take null values.

The bending moment and the axial force are distributed as follows:

[TABLE]

where $M_{\mathrm{A}},P_{\mathrm{Ax}},P_{\mathrm{Az}}$ are unknown bending moment and reaction forces in the support A. Equilibrium constraints provide additional equations for determination of the reactions at the support B. Initially, distribution of the transverse displacements will be determined from the differential equation Eq. (24):

[TABLE]

augmented with the boundary conditions:

[TABLE]

Since the differential equation is of the fourth order, only four boundary conditions can be specified at this point. This will provide a solution that involves unknown reactions in support A, $M_{\mathrm{A}}$ and $P_{\mathrm{Az}}$. To determine unknown reactions, the remaining boundary conditions in the support B, $v(L)=0,v^{\prime}(L)=-\varphi(L)=0$, have to be exploited. Determination of these can be also suitably performed by Wolfram Mathematica, but due to lengthy forms, support reactions and deflection line equations are not reported. It is just pointed out that reactions show dependence on the nonlocal parameter in the line with the discussion in Example 5.2.

As for the axial displacements, these can be determined from Eq. (22) as:

[TABLE]

together with the boundary conditions:

[TABLE]

The solution includes unknown axial reaction $P_{\mathrm{Ax}}$ at the support A. This reaction follows from $u_{0}(L)=0$.

The main results are summarized in Figs. (8-13). To make figures more easier to interpret, the length was normalized and denoted with an overbar. The minimal value for $\lambda$ serves as the limit case that corresponds to the classical Bernoulli-Euler local solution. Similarly to Example 5.2, the nanobeam behaviour tends to be more flexible until a certain threshold is reached and subsequently converge to another, more rigid curve.

The solution for the axial displacement $u(x,z)$ is graphically presented in Fig. 12 for a selected value of the nonlocal parameter. Provided distribution shows that boundary conditions are fulfilled at the left and right support. The maximal longitudinal displacement occurs at approximately 60% of the beam’s length.

Finally, Fig. 13 provides an insight into distribution of the nonlocal strain $\varepsilon(x,z)$. Compared to the more simpler case in Example 5.2, Fig. 7, the normal strain shows significantly more complex pattern. The accompanying stress field is presented in Fig. 14. Note that the local solution for stresses would provide the same distribution as for the normal strain with a difference being only due to multiplication by Young’s modulus $E$. Difference in the distribution pattern in Fig. 13 and Fig. 14 clearly indicates the nonlocal behaviour.

6 Conclusions

A new thermoelastic integral model for nanobeams has been developed by adopting the nonlocal elasticity theory by G. Romano and R. Barretta [25]. Effectiveness of the presented methodology has been tested by examining selected case-studies.

Some closing remarks are collected below.

The new stress-driven integral thermoelastic law, equipped with the averaging kernel introduced by Helmholtz, is equivalent to a set of differential equations and to suitable higher-order constitutive boundary conditions. 2. 2.

The example section also demonstrates that thermoelastic nonlocal solutions do not necessarily lead toward gradually stiffer solution for increasing values of the nonlocal parameter. Such behaviour has been displayed in statically indeterminate structures. Otherwise, solutions show monotonically increasing stiffness like in isothermal cases. 3. 3.

Since the present method is stress-driven, the dependence of the boundary conditions on the constitutive behaviour is circumvented in the statically determine problems, as demonstrated in earlier researches within isothermal framework. As discussed in Example 5.2 and further confirmed in Example 5.4, in the statically indeterminate cases the reactions still show dependence on the nonlocal parameter due to coupling of the local (thermal expansion/contraction) and nonlocal (mechanical) terms. 4. 4.

Having latter remark in mind, proposed formulation can be considered as the mixed local-nonlocal approach.

7 Acknowledgements

This work has been supported in part by Croatian Science Foundation under the project no. 6876 - Assessment of structural behaviour in limit state operating conditions and in part under the project ReLUIS - Italian Department of the Civil Protection 2017. This support is gratefully acknowledged.

References