Nonlocal integral thermoelasticity: a thermodynamic framework for functionally graded beams

TL;DR

This paper develops a thermodynamically consistent nonlocal integral model for functionally graded beams under thermal and mechanical loads, revealing complex behaviors influenced by nonlocality and temperature effects at small scales.

Contribution

It introduces a novel stress-driven integral formulation for thermoelastic analysis of functionally graded beams, incorporating thermodynamic consistency and boundary conditions.

Findings

Nonlocality and temperature significantly affect beam behavior.

Complex structural responses occur in nonisothermal indeterminate problems.

The model captures the shift of the neutral surface due to thermal effects.

Abstract

An effective nonlocal integral formulation for functionally graded Bernoulli-Euler beams in nonisothermal environment is developed. Both thermal and mechanical loadings are accounted for. The proposed model, of stress-driven integral type, is shown to be governed by a thermodynamically consistent differential problem with proper constitutive boundary conditions. The new thermoelastic strategy is illustrated by investigating a set of examples. It is demonstrated that in nonisothermal statically indeterminate problems rather complex structural behaviours can appear and that both the shift of the neutral surface and nonlocality have a dominating influence at small-scales.

Figures

Equations

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

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

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

$$w^{(1)}(x)=-\phi (x).$$

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

$$\epsilon ={\epsilon }_{\mathrm{T}}+{\epsilon }_{\mathrm{M}},$$

$${\epsilon }_{\mathrm{M}}=\epsilon -\alpha \Delta \theta =u_0^{(1)}-w^{(2)}(z-z_0)-\alpha \Delta \theta .$$

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

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

$$\begin{array}[]{l}\frac{{\sigma}}{E}=-L_{\mathrm{c}}^{2}{(\varepsilon-\alpha\Delta\theta)}^{(2)}+\varepsilon-\alpha\Delta\theta,\end{array}$$

$$\begin{array}[]{l}\sigma=E\left(-L_{\mathrm{c}}^{2}\varepsilon_{\mathrm{M}}^{(2)}+\varepsilon_{\mathrm{M}}\right),\end{array}$$

$$\begin{array}[]{l}\frac{{\sigma}}{E}=-L_{\mathrm{c}}^{2}{(u_{0}^{(3)}-w^{(4)}(z-z_{0})-(\alpha\Delta\theta)^{(2)})}\\ +u_{0}^{(1)}-w^{(2)}(z-z_{0})-\alpha\Delta\theta.\end{array}$$

$$\begin{array}[]{l}\left.L_{\mathrm{c}}\left(\varepsilon^{(1)}-(\alpha\Delta\theta)^{(1)}\right)-(\varepsilon-\alpha\Delta\theta)=0\right|_{\text{at }(0,z)},\\ \left.L_{\mathrm{c}}\left(\varepsilon^{(1)}-(\alpha\Delta\theta)^{(1)}\right)+(\varepsilon-\alpha\Delta\theta)=0\right|_{\text{at }(L,z)}.\\ \end{array}$$

$$\begin{array}[]{l}L_{\mathrm{c}}\left(u_{0}^{(2)}-w^{(3)}(z-z_{0})-(\alpha\Delta\theta)^{(1)}\right)\\ -(u_{0}^{(1)}-w^{(2)}(z-z_{0})-\alpha\Delta\theta)=\left.0\right|_{x=0},\\ L_{\mathrm{c}}\left(u_{0}^{(2)}-w^{(3)}(z-z_{0})-(\alpha\Delta\theta)^{(1)}\right)\\ +(u_{0}^{(1)}-w^{(2)}(z-z_{0})-\alpha\Delta\theta)=\left.0\right|_{x=L}.\\ \end{array}$$

$$u_0^{(1)}{\int }_{\Omega }E\mathrm{d}A-w^{(2)}{\int }_{\Omega }E(z-z_0)\mathrm{d}A-{\int }_{\Omega }E\alpha \Delta \theta \mathrm{d}A={\int }_0^L{\varphi }_{\mathrm{c}}(x-\xi )N(\xi ,z)\mathrm{d}\xi ,$$

$$\begin{array}[]{l}L_{\mathrm{c}}\left(u_{0}^{(2)}\int_{\Omega}E\mathrm{d}A-w^{(3)}\int_{\Omega}E(z-z_{0})\mathrm{d}A-\int_{\Omega}E(\alpha\Delta\theta)^{(1)}\mathrm{d}A\right)\\ -(u_{0}^{(1)}\int_{\Omega}E\mathrm{d}A-w^{(2)}\int_{\Omega}E(z-z_{0})\mathrm{d}A-\int_{\Omega}E\alpha\Delta\theta)\mathrm{d}A=\left.0\right|_{x=0},\\ L_{\mathrm{c}}\left(u_{0}^{(2)}\int_{\Omega}E\mathrm{d}A-w^{(3)}\int_{\Omega}E(z-z_{0})\mathrm{d}A-\int_{\Omega}E(\alpha\Delta\theta)^{(1)}\mathrm{d}A\right)\\ +(u_{0}^{(1)}\int_{\Omega}E\mathrm{d}A-w^{(2)}\int_{\Omega}E(z-z_{0})\mathrm{d}A-\int_{\Omega}E\alpha\Delta\theta)\mathrm{d}A=\left.0\right|_{x=L}.\\ \end{array}$$

$$u_0^{(1)}{\int }_{\Omega }E(z-z_0)\mathrm{d}A-w^{(2)}{\int }_{\Omega }E(z-z_0{)}^2\mathrm{d}A-{\int }_{\Omega }E\alpha \Delta \theta (z-z_0)\mathrm{d}A={\int }_0^L{\varphi }_{\mathrm{c}}(x-\xi )M(\xi ,z)\mathrm{d}\xi ,$$

$$\begin{array}[]{l}L_{\mathrm{c}}\left(u_{0}^{(2)}\int_{\Omega}E(z-z_{0})\mathrm{d}A-w^{(3)}\int_{\Omega}E(z-z_{0})^{2}\mathrm{d}A-\int_{\Omega}E(\alpha\Delta\theta)^{(1)}(z-z_{0})\mathrm{d}A\right)\\ -(u_{0}^{(1)}\int_{\Omega}E(z-z_{0})\mathrm{d}A-w^{(2)}\int_{\Omega}E(z-z_{0})^{2}\mathrm{d}A-\int_{\Omega}E\alpha\Delta\theta(z-z_{0})\mathrm{d}A)=\left.0\right|_{x=0},\\ L_{\mathrm{c}}\left(u_{0}^{(2)}\int_{\Omega}E(z-z_{0})\mathrm{d}A-w^{(3)}\int_{\Omega}E(z-z_{0})^{2}\mathrm{d}A-\int_{\Omega}E(\alpha\Delta\theta)^{(1)}(z-z_{0})\mathrm{d}A\right)\\ +(u_{0}^{(1)}\int_{\Omega}E(z-z_{0})\mathrm{d}A-w^{(2)}\int_{\Omega}E(z-z_{0})^{2}\mathrm{d}A-\int_{\Omega}E\alpha\Delta\theta(z-z_{0})\mathrm{d}A)=\left.0\right|_{x=L}.\\ \end{array}$$

$${\sigma }^{(1)}=0.$$

$$\rho \dot{e}+q^{(1)}=\sigma +\rho r+\mathcal{P},$$

$$\rho \dot{\eta }\ge \frac{\rho r}{\theta }-{\left(\frac{q}{\theta }\right)}^{(1)}.$$

$$g(\sigma ,\theta )=e-\eta \theta -\sigma \epsilon .$$

$$\psi =g+\sigma \epsilon .$$

$$-\rho \dot{g}-\dot{\sigma }\epsilon -\rho \eta \dot{\theta }-\frac{q}{\theta }{\theta }^{(1)}\ge 0.$$

$$\dot{g}={\partial }_{\sigma }g\dot{\sigma }+{\partial }_{\theta }g\dot{\theta },$$

$$-(\epsilon +{\partial }_{\sigma }g)\dot{\sigma }-\rho ({\partial }_{\theta }g+\eta )\dot{\theta }-\frac{q}{\theta }{\theta }^{(1)}\ge 0.$$

$$\epsilon =-{\partial }_{\sigma }g.$$

$$\eta =-{\partial }_{\theta }g.$$

$$-\frac{q}{\theta }{\theta }^{(1)}\ge 0.$$

$$\begin{array}[]{c}\rho g(\sigma,\theta)=-\frac{1}{2E}\sigma^{2}-\alpha\Delta\theta\sigma.\end{array}$$

Full text

††journal: Composite Structures

1 Introduction

Beam models that can capture small size behaviour, typically observed in very small devices, is a very active research field nowadays. As small beams, we presume structures having dimensions of micrometer or nanometer order. Size effects (i.e. nonlocal phenomena) are observed since the forces irrelevant at the classical engineering scale (where structures are sized in centimetres or meters) are now important in such small structures. Increased technological developments at microscale and nanoscale require more accurate mechanical solutions, fuelling the increased interest of the research community [1, 2, 3].

Initial contributions to the solution of these problems were based on the gradient formulation. Although some paradoxical results were reported rather early in literature [4, 5], such an approach remains a popular choice. The paradox has puzzled scientists for quite some time, and some attempts based on the strain-driven integral formulations were proposed in literature. Unfortunately, some strain-driven models suffer from the problem of fulfilling equilibrium conditions [6]. On the contrary, integral formulations of the stress-driven type [7] are able to properly address issues relating to the fulfilment of equilbrium conditions and paradoxes debated in the scientific community [8, 9]. Moreover, it is pointed out that, at least for a certain class of kernel function, the stress-driven integral formulation can be made equivalent to differential problems [10, 11, 12]. To ensure this equivalence, suitable constitutive boundary conditions must be set in the correct manner. These are different in isothermal and nonisothermal problems.

Nonisothermal environmental conditions can have a profound influence on the mechanical response of nonlocal beams. In the vibrational context, temperature has a significant influence on the frequency shift [13]. Thus, a good part of research efforts on nonlocal beams in nonisothermal environments are devoted to dynamical analysis [14]. Composite beams in this context are analysed in [15, 16, 17, 18, 19]. These analyses tend to account for other phenomena and become rather complex, like in the thermo-electro-mechanical analysis of functionally graded size dependent beams [20]. Non-probabilistic uncertainty modelling for vibration and buckling of the FG nanobeams in nonisothermal conditions is considered in [21]. In the statical analysis, buckling of nonlocal functionally graded beams under thermal loading is frequently addressed, starting from the early contribution in [22] till more recent contributions [23, 24, 25, 26, 27]. Mainly nonlocal influence on the mechanical part of the problem is investigated, but buckling caused by the size effect on heat conduction is also analysed [28]. Another recent contribution [29] points out the difference in original and simplified boundary conditions in vibrational and buckling analyses of FG beams in nonisothermal environments. Original boundary conditions do include the thermal moment at the flexural boundary conditions. The obtained results clearly demonstrate significant differences of two types of boundary conditions. Bending of thermoelastic nonlocal beams were analysed in [30, 31]. Majority of these results are based on gradient formulations, both for Bernoulli-Euler and Timoshenko beams. Only a few solutions based on the integral formulation are available, see [32] for an example.

It is noteworthy that in composite beams the neutral surface does not coincide with the geometrically neutral surface [33, 34]. This is a direct consequence of the functional grading of beam’s material. Such effects does bring additional level of complexity to the modelling procedure. Even in the classical, local case the most frequent approach is to disregard these effects. Nevertheless, some results including this shift in local problems do exists. For the present research, nonisothermal model [29] is of interest. But, when the nonlocal procedures are concerned, inclusion of the shift of neutral surface is much less common; for the isothermal strain gradient formulation see [35, 36, 37, 38] and for the integral based [39]. For the gradient based thermoelastic beams including hygro effects, see [15]. In the case of stress-driven nonlocal thermoelasticity, these are practically non-existent. This is surprising, since its level of influence on the results should be at least the same as the small size effects. Thus, this paper aims to include such effects into a thermoelastic formulation.

To summarize, novelty of the present paper is inclusion of the neutral surface shift into the formulation by extending previous results obtained for the nonlocal nonisothermal Bernoulli-Euler formulation [32] in the case of homogeneous beam. As a consequence, the formulation obtained in such manner is suitable for beams made of functionally graded materials. It is also shown that the strict thermodynamic framework is imperative in these derivation. Starting from the first and second law of thermodynamics, a suitable Gibbs potential is developed including all constraints arising from the stress-driven integral nature of the model. Minimization of this newly proposed potential yields the system of ordinary differential equations governing the mechanical process. The proposed formulation is thoroughly tested on four examples.

The paper is organized as follows: in order to provide the notation and a stage for further developments, the Bernoulli-Euler beam kinematics is presented at the start. The next section describes the transformation from the integral law to the equivalent differential formulation. Subsequently, to accommodate this formulation, a strict thermodynamic framework is introduced. A cornerstone of the formulation is the new potential; minimization of the Gibbs potential provides governing equations as well as the complete set of boundary conditions. The example section is followed by the conclusions that summarize main findings and close the paper.

2 Beam kinematics in the nonisothermal setting

To keep elaborations as simple as possible, plane bending of an initially straight Bernoulli-Euler beam will be considered. The longitudinal axis of the beam is denoted as $x$, while the bending is assumed to take place in the $x-z$ plane. The z-axis is assumed to originate at the geometric centroid of the cross-section, so that the first moment of area vanishes $\int_{\Omega}z\mathrm{d}A=0$. The beam is made of functionally graded (FG) material, where grading does take place along $z$ axis.

In the case of a FG beam, the physical neutral surface, in which the normal stress vanishes, does not coincide with the geometrical middle surface. The issue is recently raised in [29] for thermally loaded FG local beams. Thus, in calculation of the axial displacement, the shift of neutral surface from the geometrical middle surface $z_{0}$ must be accounted for. It will be demonstrated that $z_{0}$ can vary along the beam if $z_{0}=z_{0}(x)$. Nevertheless, the present research will not address such problems.

The central task is to determine the Cartesian components of the displacement vector field. Having in mind the remark about the shift of neutral surface, these are [40]:

[TABLE]

Above, $\varphi(x)$ denotes the angle of rotation of the cross section. The longitudinal displacement $u_{x}$ is composed of two parts. The first part $u_{0}$ represents the average displacement of the cross section, defined as the integral of the displacement over the cross sectional domain $\Omega$:

[TABLE]

where $A$ denotes the cross sectional area. In the Bernoulli-Euler formulation cross sections remain orthogonal to the beam axis. This implies vanishing shear strains $\gamma_{xz}$ as the sum of the derivative of longitudinal and transversal displacement:

[TABLE]

providing the link between the transversal displacement and angle of rotation:

[TABLE]

The apex (n) denotes the $n$-th derivative with respect to the longitudinal coordinate $x$. The axial strain in a point follows from the differentiation of the axial displacement field with respect to $x$:

[TABLE]

Such a kinematic framework accounts for both isothermal and nonisothermal effects.

It remains to explicitly introduce the temperature effects. In order to do so, the normal strain is additively separated into a thermal part $\varepsilon_{\mathrm{T}}=\alpha\Delta\theta$ and a mechanical part $\varepsilon_{\mathrm{M}}$:

[TABLE]

where $\alpha$ is the coefficient of thermal expansion and $\Delta\theta(x,z)$ is the temperature change field. With Eqs. (5,6), the mechanical strain is:

[TABLE]

In the local formulation, the mechanical strain is the source of normal stress, whereas in the nonlocal formulation its determination is not that much straightforward.

3 Thermodynamic preliminaries for functionally graded beams

3.1 Motivation

The cornerstone of the FG formulation is the stress-driven nonlocality defining constitutive law for the mechanical part of axial strain:

[TABLE]

where the kernel function $\phi_{\mathrm{c}}({x)}$ is

[TABLE]

${\sigma}$ is the axial stress and $E$ is the Euler-Young’s modulus. The characteristic length is defined as $L_{\mathrm{c}}=cL$, i.e. as the product of the small-size parameter $c$ and the beam’s length $L$. Note that such assumption accounts only for nonlocality with respect to the longitudinal coordinate. 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 that are situated at a larger distance.

It should be emphasized that the precise value of the small size parameter still remains an open question. Preferably, the value should be determined experimentally, but only rare attempts to quantify size-dependent behaviour can be found in the literature, see [41, 42, 43] for an insight. Available numerical calculations [44, 45, 46, 47, 48, 49] also indicate that further work remains to be done in order to solve the issue.

If one wishes to obtain stresses from known strains, it can shown that the solution and necessary boundary conditions follow from the differential formulation [50, 32]:

[TABLE]

or alternatively in a more compact form:

[TABLE]

cf. Example 5.2 [51] in the case of gradient based methods. In the local case, the latter equation collapses to the usual form $\frac{{\sigma}}{E}=\varepsilon_{\mathrm{M}}=\varepsilon-\alpha\Delta\theta$. Application of Eq. (7) that defines the link between strain and displacement, transforms the stress into:

[TABLE]

In contrast to the widely employed gradient formulations, here additional constraints in the form of boundary conditions have to be imposed at $x=0$ and $x=L$ in order to assure solution of the constitutive integral model [50]:

[TABLE]

These constraints assume existence and uniqueness of the solution Eq. (10) [32]. Obviously, constraints relate kinematical quantities and imply that standard kinematical boundary conditions cannot be chosen arbitrarily. These are known as the constitutive boundary conditions and were introduced in the nonlocal beam formulation in [6, 7]. Consequently, the stress-driven integral formulation can be made equivalent to the gradient one if the boundary conditions Eq. (13) are respected.

Alternatively, by the application of displacement and rotation, Eq. (5), the above constraints become:

[TABLE]

The stress-driven model valid in a beam’s point should be now homogenized for an arbitrary cross section. Hence, Eq. (8) is multiplied by the Euler-Young’s modulus (note that it must be $E\neq E(x)$ in order to do so), integrated over the cross section and finally the equilibrium equation $\int_{\Omega}\sigma\mathrm{d}A=N$ is applied. Introducing displacements Eq. (5), this provides:

[TABLE]

while the accompanying constitutive boundary conditions become:

[TABLE]

The same reasoning as for the equilibrium equation in the axial direction can be now applied to bending, $\int_{\Omega}\sigma(z-z_{0})\mathrm{d}A=M$. Multiplication of Eq. (8) by $(z-z_{0})$ yields:

[TABLE]

and

[TABLE]

This constitutive model must be now placed in the strict thermodynamic framework.

3.2 Continuum mechanics fundamentals of beam-like structures

Having introduced the cornerstone of the stress-driven nonlocal model, the proper thermodynamic framework specialized to beams is provided. For this purpose, a brief recapitulation of balance equations from continuum mechanics is presented. Due to one dimensional nature of the problem, vectorial and tensor quantities collapse to the scalar form. The balance of momentum in a Bernoulli-Euler beam where inertial and volume forces are deemed to be insignificant can be written as:

[TABLE]

Likewise, the first law of thermodynamics is [52]:

[TABLE]

while the second law is:

[TABLE]

Above, $\rho(x,t)$ is the specific mass, $e(x,t)$ is the internal energy, $q(x)$ is the heat flux, $d$ is the rate of deformation, $r(x,t)$ is the heat source per unit mass, $\theta(x,t)$ is temperature, $\mathcal{P}$ is the rate of energy supply from the neighbourhood and $\eta(x,t)$ is the specific entropy per unit mass. The origin of the rate of energy supply $\mathcal{P}$ is nonlocality. Note that for the whole beam it must be $\int_{\mathcal{B}}\mathcal{P}\mathrm{d}V=0$ in order to fulfil the balance of energy, see [51] for more details. As usual, a superimposed dot represents differentiation with respect to time. We now postulate existence of the Gibbs potential $g(\sigma,\theta)$ [53]:

[TABLE]

This introduces the link between the Helmholtz free energy and Gibbs potential as:

[TABLE]

Although generally valid, the latter equation does not imply that an equivalent Helmholtz function can be obtained from the Gibbs functions in all cases (and vice-versa). As demonstrated in [54], such transformation is not always possible. Thus, for some type of problems, the Helmholtz free energy function is more appropriate while for other the Gibbs function is the preferred choice.

With the Gibbs potential at hand, it can be shown that the second law of thermodynamics takes the form:

[TABLE]

Introducing the rate of the Gibbs potential as

[TABLE]

it is obtained:

[TABLE]

Now, standard arguments are invoked [55]. Since the last equation has to hold for all arbitrary processes, it is obtained:

[TABLE]

Finally, with Eq. (27) at hand, consider a thermoelastic process with homogeneous temperature field. Inequality (26) yields the definition of entropy:

[TABLE]

In purely thermal processes in which deformation does not take place it is:

[TABLE]

This effectively puts constraints on the allowable direction of the heat flow. Since the present research will not address problems involving heat conduction, the last equation will not be used in the forthcoming elaborations.

3.3 Prototype model

Having set up the general thermodynamic framework for beams, some constitutive assumptions must be now introduced. In particular, the constitutive behaviour by means of the Gibbs function and a functional for determination of kinematical quantities is provided.

3.3.1 Gibbs potential

In order to put the proposed nonlocal stress-driven formulation within the thermodynamic framework described above, the Gibbs potential $g$ of the following form is used in the local formulation [56]:

[TABLE]

In this particular case the above function should be modified in order to accommodate the nonlocal character of the model Eq. (8):

[TABLE]

Above, the Euler-Young’s modulus is assumed to be a function of the transverse coordinate, $E=E(z)$. It can be easily verified that the above choice of the Gibbs function reflects Eqs. (8) if introduced into Eq. (27).

3.3.2 Governing potential for beams

To provide a variationally consistent framework, the existence of the governing potential is postulated:

[TABLE]

The internal potential $U$ is defined for the complete domain of the beam $\mathcal{B}$ for a reversible steady state thermal process ($\dot{\theta}=0$) as:

[TABLE]

where Eq. (27) was applied. The potential of external forces $W$ assumes the standard form:

[TABLE]

where $\mathcal{N}$ are the external axial forces while $\mathcal{M}$ are bending moments at the beam end’s. Note that constitutive boundary conditions Eqs. (16,18) arising from integrability constraints are indirectly included.

Finally, the solution follows from the minimum total potential energy principle:

[TABLE]

Required stationary conditions are discussed in derivations that follow. Finally, the variational problem at hand can be alternatively solved by the application of the minimum total complementary energy principle.

3.3.3 Shift of the neutral surface

Yet remains to define the shift of the neutral surface from the centroid of the cross section. In the case of homogeneous beams, the beam’s cross section rotates around the $y$ axis that is perpendicular to the plane of deformation. This axis is situated at the centroid of the cross section. It is understood that the lines parallel to the longitudinal axis laying in the neutral surface are bent, but their length does not change what results in vanishing normal stresses. However, if the beam is functionally graded, the situation is changed. Variation of beam stiffnesses can cause the shift of the neutral surface $z_{0}$. The section at hand will demonstrate how to determine this shift.

In the local thermoelastic problem, the displacement of any point of the beam can be considered to consists of two parts: translational and rotational part. Likewise, the axial strain of mechanical origin in a longitudinal cross section defined by the coordinate $z$ is composed of three parts, Eq. (7). Translation is related to the uniform elongation term $u_{0}^{(1)}$ and thermal elongation $\alpha\Delta\theta$, while the rotation due to bending causes $-w^{(2)}(z-z_{0})$. In the neutral surface, axial displacement due to rotation (and consequently strain due to bending) must vanish. This can be used to determine the shift $z_{0}$. Now, the stresses in an arbitrary cross section $x$ are calculated as:

[TABLE]

Since in the neutral surface stresses (and strains) due to rotation of the cross section should vanish it is obtained:

[TABLE]

Now, integration of the above equation over the cross section provides the shift $z_{0}$:

[TABLE]

Note that in the case of $E=E(x,z)$, the above shift becomes a function of the axial coordinate, $z_{0}=z_{0}(x)$ significantly complicating the formulation. However, in some special cases simpler solutions can be obtained. For instance, let the Euler-Young’s modulus can be expressed as $E(x,z)=E_{x}(x)E_{z}(z)$. The advantage of such assumption is obvious when $E$ is introduced into Eq. (38):

[TABLE]

conveniently transforming $z_{0}$ to a constant. Obviously, the identical equation is obtained for the case $E=E(z)=E_{z}$.

If a beam is homogeneous ($E=\text{const.}$) or $E=E(x)=E_{x}$, the shift of the neutral surface $z_{0}$ is zero, since Eq. (38) collapses to

[TABLE]

Note that the first moment of area still vanishes in the FG formulation, $\int_{\Omega}z\mathrm{d}A=0$.

If the nonlocal formulation, the starting counterpart is Eq. (8). Suppose again that $E(x,z)=E_{x}(x)E_{z}(z)$. This gives:

[TABLE]

Following the same reasoning as in the local formulation, at the neutral surface strain contribution from the rotation term must vanish, again providing:

[TABLE]

repeating the result from the local case.

3.4 Stationarity conditions

3.4.1 Stationarity with respect to $\delta_{{u}_{0}}\Pi$

This section deals with the stationary conditions that provide governing equations and boundary conditions. First term $\delta_{{u}_{0}}\Pi=0$ is obtained as follows. Starting from the total energy potential Eq. (32) and strain Eq. (5), it is obtained:

[TABLE]

Integration by parts with respect to the longitudinal coordinate $x$ provides:

[TABLE]

The first term will provide boundary conditions as

[TABLE]

It is emphasized that constitutive boundary conditions Eq. (16) accompany above listed standard boundary conditions.

The second term in Eq. (44) gives in the same manner:

[TABLE]

or having in mind arbitrariness of the virtual axial displacement:

[TABLE]

Due to Eq. (11), $\sigma^{(1)}$ is:

[TABLE]

Explicit forms for $\varepsilon_{\mathrm{M}}^{(1)}$ and $\varepsilon_{\mathrm{M}}^{(3)}$ will be introduced later. Alternatively, Eq. (47) can be also rewritten as:

[TABLE]

if the equilibrium equation in the axial direction is employed.

3.4.2 Stationarity with respect to $\delta_{{w}}\Pi$

Second stationary condition $\delta_{{w}}\Pi=0$ gives:

[TABLE]

The same steps as in the first stationary condition are followed. Integrating by parts twice yields:

[TABLE]

The first two terms with Eq. (12) again provide boundary conditions. The first boundary condition exploits arbitrariness of $\delta w^{(1)}$ at $x\in\{0,L\}$:

[TABLE]

and the second one $\delta w$ at $x\in\{0,L\}$:

[TABLE]

The constitutive boundary conditions Eq. (18) arising from the integrability conditions must be included as well.

The third and fourth term in Eq. (51) gives:

[TABLE]

or

[TABLE]

The first derivative of stress is defined by Eq. (48), while the second one is:

[TABLE]

Like in the axial case, Eq. (55) can be transformed by means of the equilibrium equation into:

[TABLE]

4 Beam displacements of nonlocal FG beams

Equations defining beam displacements are derived for functional grading in one direction of the Euler-Young’s modulus $E=E(z)$. The coefficient of thermal expansion and temperature are of a general form: $\alpha=\alpha(x,z)$, $\theta=\theta(x,z)$. At the same time, following notation is introduced:

[TABLE]

In the present case, the shift of the neutral surface is not zero. The first and second derivative of the stress $\sigma^{(1)},\sigma^{(2)}$ now have nonlocal character, Eq. (10):

[TABLE]

With above result at hand, the governing equation for the axial displacement follows from Eq. (47) as

[TABLE]

where some terms were cancelled due to the property of the neutral surface. With the notation Eq. (58), the equation describing axial deformation is:

[TABLE]

Boundary conditions again follows from Eq. (45) as:

[TABLE]

and from Eq. (16):

[TABLE]

The governing equation for bending is derived from Eq. (55):

[TABLE]

Introduction of the second derivative of stress, Eq. (59)2 and accounting for the neutral surface effects, it is obtained:

[TABLE]

cf. Remark 7.1 in [7] for the isothermal case. Thus, in the contrast to the standard fourth order differential equation, the present sixth order case includes constitutive boundary conditions Eq. (68) into the formulation.

The boundary conditions are defined by Eqs. (52, 53):

[TABLE]

and the second one:

[TABLE]

along with conditions Eq. (18)

[TABLE]

With above differential equations and accompanying boundary conditions the formulation is completely defined.

5 Examples

5.1 Homogeneous doubly clamped bar heated by the uniform thermal field

The introductory example deals with a nonlocal doubly clamped bar of length $L$ subjected to the uniform temperature field $\Delta\theta$. The similar example was analysed in [32] (Example 5.2), but with formulation that is derived slightly differently than in the present example. In particular, the approach chosen in the latter paper uses differential equation of third order for describing axial displacements, while here fourth order differential equation is employed. Naturally, this requires different number of boundary conditions. The bar is not functionally graded and material properties are defined by constant $E,\alpha$. The cross section is rectangular with dimensions $b$ and $h$.

Thus, the problem that should be solved is governed by Eq. (61):

[TABLE]

The constitutive boundary conditions are:

[TABLE]

At the left end, $u_{0}(0)=0$ was prescribed. The boundary condition at the right end follows from Eq. (62):

[TABLE]

where $P$ denotes the unknown axial reaction in the support. The solution of the problem coincides with the one in [32]:

[TABLE]

The unknown support reaction follows from the remaining boundary condition $u_{0}(L)=0$ as:

[TABLE]

or after straightforward algebraic manipulation:

[TABLE]

This confirms the equivalence of both approaches.

5.2 Axial thermal deformation of a homogeneous doubly clamped bar heated by the longitudinally varying temperature field

The second example deals with a nonlocal doubly clamped bar of length $2L$ subjected to the varying temperature field $\Delta\theta=-\theta_{0}\left(1-\frac{x^{2}}{L^{2}}\right)$. The origin of the coordinate system is positioned at the midpoint of bar, so the endpoints are situated at $x=-L$ and $x=L$. Thus, $L_{c}=2L\;c$. The present example is chosen in order to compare the solution of the strain gradient model presented in [51], Example 5.2 to the integral stress-driven model introduced in this manuscript. Material properties and cross section dimensions are similar like in the previous example.

Governing differential equation is again defined by Eq. (61):

[TABLE]

where $\Phi^{\mathrm{N}}_{3}=0$, $\Phi^{\mathrm{N}}_{1}=2\alpha E\theta_{0}\frac{A}{L^{2}}x$ and $q_{x}=0$ were used. The governing equation is the same for the stress-driven formulation and the gradient one. The general solution of the equation is:

[TABLE]

where $C_{i},i=1,2,3,4$ are integration constants. To obtain these, boundary conditions must be enforced. In both formulations, displacements at the end of the bar must be zero:

[TABLE]

The rest of boundary conditions are different in each formulation. In the gradient one these are:

[TABLE]

In the stress-driven integral formulation, the boundary conditions that must be fulfilled follow from Eq. (16)

[TABLE]

The gradient formulation provides axial displacement as:

[TABLE]

while the integral one gives:

[TABLE]

Support reactions are determined as $\left.P=\int_{\Omega}\sigma\mathrm{d}A\right|_{x=L}$:

[TABLE]

All results are graphically presented for $L=1,\alpha=0.1,\theta_{0}=10,c=0.2$, $b=1$, $h=1$. Dependence of support reactions on the small size parameter is presented in the logarithmic plot in Fig. 1. Obviously, the gradient approach is governed by the quadratic function and results in much higher reaction forces than the integral one as the small size parameter is increased. The integral formulation results in a more complex dependence, most easily graphically interpreted. This formulation shows much slower increase than the gradient method reaction.

Solutions for axial displacements also differ. Distribution of the axial displacement $u_{0}$ shows almost identical solutions for the smallest values of $c$, Fig. 2. As the small size parameter is increased, solutions obtained with the gradient method show monotonic increase in beam’s stiffness. Displacements at the beam ends are obviously identical due to boundary conditions Eq. (77) for both formulations. Unlike the gradient solution, the integral approach initially exhibits stiffening, but after a certain threshold the displacements again increase. This can be explained by the different last terms in parentheses of Eqs. (80, 81). Such behaviour was already noticed earlier, see [32] for further clarifications. The strain distribution, Fig. 3 also manifests differences between two methods. The origin of such behaviour is the difference in boundary conditions Eqs. (78, 79). Since the gradient approach enforces the second derivative of the axial displacement to be zero at the beam ends (Fig. 4), this has the consequence that the tangent to curve describing the strain distribution must be horizontal at ends. This is not the case in the integral formulation.

5.3 Thermally loaded FG cantilever beam

To demonstrate dependency of results on the neutral surface shift, a cantilever FG beam is considered. Temperature field $\Delta\theta=\theta_{1}z^{2}$ is imposed on the beam. At $x=L$ the beam is loaded by the axial force $P$ and couple $M$. The coefficient of thermal expansion is $\alpha$, the Euler-Young’s modulus is linearly distributed across the height of the beam $E(z)=E_{0}+E_{1}(z+h/2)$, while the beam rectangular cross section is defined by the height $h$ and the width $b$ as usual.

As a consequence of the functional grading, the neutral surface is shifted by (Eq. (42)):

[TABLE]

The problem is defined by the set of differential equations Eqs. (61, 65):

[TABLE]

where terms $\Phi^{\mathrm{N}}_{1}$, $\Phi^{\mathrm{N}}_{3}$, $\Phi^{\mathrm{M}}_{2}$, $\Phi^{\mathrm{M}}_{4}$ vanish since the temperature does not depend on the longitudinal coordinate $x$. The set of boundary conditions follows from Eqs. (62,63,66,67,68):

[TABLE]

where $\Phi^{\mathrm{N}}_{2}$, $\Phi^{\mathrm{M}}_{3}$, $\Phi^{\mathrm{M}}_{1}$ again disappear.

Solutions of the above problem are:

[TABLE]

Note that in the absence of the external loads, $P=0$ and $M=0$, the nonlocal character of solutions disappears.

To graphically illustrate solutions, the following values are selected: $L=1,\alpha=1,\theta_{1}=0.1,b=1,h=1,P=1,M=1$. Both solutions were evaluated for the small size parameters $c\in\left\{0.0001,0.01,0.02,0.05,0.2,0.5,1.0,2.0,10\right\}$. The results are given in Figs. 5, 6. Both displacements show the stiffening effect with the increase in the small size parameter.

For the purpose of comparison, the solution for the transverse displacement for the case when the shift of the neutral surface is neglected, $z_{0}=0$ is also reported:

[TABLE]

so the transverse displacement assumes simpler form if $z_{0}=0$. The axial displacement is not affected by the neutral surface shift.

The significance of the neutral surface shift is illustrated in Fig. 7. Increase in beam’s stiffness with the increase in $\left|z_{0}\right|$ is clearly noticeable.

5.4 Thermomechanically loaded simply supported FG beam

A simply supported FG beam is loaded by the distributed axial loading $q_{x}$ and distributed transverse loading $q_{z}$. Additionally, beam is exposed to the bilinear temperature field $\Delta\theta=\theta_{0}xz$. The beam material is linearly functionally graded like in the previous example $E(z)=E_{0}+E_{1}(z+h/2)$, with the same geometry. The coefficient of heat expansion $\alpha$ is taken to be a constant.

The problem is defined by the differential equations Eqs. (61, 65)

[TABLE]

The set of boundary conditions follows from Eqs. (62,63,66,67,68):

[TABLE]

Solving the system of ordinary differential equations just described provides the axial and transverse displacement:

[TABLE]

Diagrams representing the deformation process are prepared for values: $L=1,\alpha=1,\theta_{0}=0.1,b=1,h=1,E_{0}=0.1,E_{1}=1$. Since the same cross section and material like in Example 5.3 are used, the shift of the neutral surface Eq. (83) remains the same: $z_{0}=\frac{2E_{1}I_{y}}{2E_{0}A-bh^{2}E_{1}}$. Dependence of the neutral surface shift of Euler-Young’s coefficients $E_{0},E_{1}$ is depicted in Fig. 8. Obviously, for $E_{1}=0$ the shift is $z_{0}=0$, while for larger values of $E_{0},E_{1}$ converges toward $z_{0}\approx 0.6$.

Obtained displacement distributions Eq. (90)1 show that the axial displacement is independent of the neutral shift. This is not the case with the transverse displacement, Fig. 9. As the shift of the neutral surface increases, the beam becomes more stiffer, asymptotically converging to the rigid beam.

Like before, solutions were evaluated for the same set of small size parameters $c\in\left\{0.0001,0.01,0.02,0.05,0.2,0.5,1.0,2.0,10\right\}$. The results are given in Figs. 10, 11 repeating behaviour observed in Example 5.3 - displacements demonstrate the stiffening effect as $c$ is increased. Hence, such a beam shows the standard nonlocal behaviour.

6 Conclusions

The paper at hand extended the existing nonlocal Bernoulli-Euler homogeneous beam formulation to functionally graded materials. The presented formulation is based on the consistent thermodynamic formulation emanating from a suitably selected potential. Most significant findings and results are summarized below.

A thermodynamically consistent formulation is employed to develop the proposed formulation, based on a proper Gibbs potential, the stress-driven integral model for functionally graded nanobeams. Minimization of the total energy potential provides underlying differential equations and boundary conditions governing the relevant thermoelastostatic problem.

• 2.

The shift of neutral surface is included into the formulation. It is demonstrated that this effect can have a profound influence on deformation, especially in the presence of small-size phenomena, and should not be disregarded.

• 3.

It is noticed that in the case of doubly-clamped nonisothermal beams, for the lower values of the nonlocal parameters, a softening behaviour compared to the local case is exhibited. However, as the nonlocal parameter is further increased, the beam starts to stiffen, asymptotically converging to the local solution. The similar conclusion is reached in [32] in the special framework of elastically homogeneous materials.

7 Acknowledgements

This work has been fully supported by the University of Rijeka under the project number uniri-tehnic-18-37 - Mechanical behaviour of nanostructures. This support is gratefully acknowledged

8 Data availability

Apart from the analytical results presented in the manuscript, there are none other data appropriate for sharing.