1 Introduction

When the internal characteristic lengths are comparable to the dimensions of structures, such as nanotubes and nanobeams, size effects manifest themselves and become significant^[1]. Such problems can be solved experimentally^[2] or by theoretical molecular/atomic modeling^[3-4]. Since both methods are quite expensive, time consuming, and requiring substantial computational resources, an alternative atomistic-continuum modeling is needed and might be of great help.

Nonlocal elasticity provides a potential way of describing the connection between microstructures and mechanical properties of materials, which can date back to Kröner^[5], Krumhansl^[6], and Kunin^[7], and was further developed by Eringen and Edelen^[8] and Eringen^[9-11] in the framework of thermodynamics. Differing from classical local continuum mechanics, the nonlocal theory abandons the classical localized hypothesis and introduces nonlocal residuals in the equilibrium equations^[11]. Besides, the long-range interactions and size effects are considered by introducing a nonlocal kernel function in the constitutive equation.

The constitutive equations commonly used are Eringen's nonlocal integral equation and differential equation, of which the latter is derived from the former for a specific class of kernel functions. Due to tremendous difficulties in solving integral equations, the nonlocal differential form is widely used.

Peddieson et al.^[12] first applied Eringen's nonlocal differential model to nanotechnology, and studied the case of a cantilever beam subject to distributed or concentrated loads, predicting that devices in nano-scales, rather than micro-scales, would exhibit nonlocal effects (or size effects). Especially, it is revealed that no size effect would happen for a nano-cantilever beam subject to a concentrated force. Challamel and Wang^[13] noticed this phenomenon, and introduced a modified lower-order nonlocal model (or gradient elasticity model) trying to resolve this paradox. Using this modified model, the research showed that when the regularization parameter is larger than 1, an increase in the small length scale parameter (or size factor) would reduce the deflection; while the deflection would increase if the regularization parameter is smaller than 1. Li^[14] studied the static and dynamic properties of nanobeams using Eringen's differential constitutive equation and obtained a softening effect, that is, an increase in the size factor would induce a growth in deflection but a decline in the fundamental frequency. However, the results from the strain gradient theory and experiments given by Lam et al.^[15] showed that the stiffness of the beam increases when the ratio of beam thickness to length scale parameter decreases. Furthermore, Li and Chou's simulation^[16] showed that the values of fundamental frequency of single- and double-walled carbon nanotubes calculated from continuum modeling are 40%~60% lower than those from atomistic modeling, which implies a stiffening effect too.

To resolve this conflict, Lim^[17] provided an effective nonlocal stress model, where the high-order equilibrium equation and high-order boundary conditions were deduced. Then, the solution of an Euler-Bernoulli beam subject to a tip load displays both size effect and a stiffening effect, that is, the deflection decreases as the size factor increases.

Li et al.^[18] started from the two-dimensional nonlocal differential constitutive equation for the plane-stress state and then derived the governing equation of nanobeams under certain assumptions. Solving a nano-cantilever beam problem, it is found that no nonlocal size effect shows up when the beam is subject to a concentrated force alone. In addition, the equivalent stiffness of a nanostructure predicted by the nonlocal theory may be larger or smaller than that given by the classical theory, depending on the types of applied loads.

It should be noticed that the above nonlocal works are all based on Eringen's differential constitutive equation, and such reported inconsistencies might indicate that one needs to study the problem based on the integral constitutive equation.

As to Eringen's integral constitutive equation, the pure nonlocal form and later proposed two-phase local/nonlocal form were commonly used. Polizzotto^[19] extended three variational principles of classical elasticity, that is, the total potential energy, the complementary energy, and the Hu-Washizu principles, to Eringen's two different nonlocal integral models, and constructed a standard finite element method (FEM) framework. Khodabakhshi and Reddy^[20] presented a general finite element formulation for the two-phase local/nonlocal elastic model, and solved the boundary-value problems of an Euler-Bernoulli beam with different types of boundary conditions and loads, showing a softening effect.

Fernández-Sáez et al.^[21] studied the inconsistencies between Eringen's original nonlocal integral model and differential model, which were once supposed to be equivalent. They found that additional conditions need to be satisfied when converting the integral equation to its differential form. Otherwise, the two equations would not be equivalent mathematically. Based on the work of Fernández-Sáez et al.^[21], Tuna and Kirca^[22] derived the closed-form analytical solutions for static bending of Euler-Bernoulli and Timoshenko beams by transforming Fredholm integral equation of the first kind to Volterra integral equations and applying the Laplace transformation. However, their work received a disagreement from Romano and Barretta^[23]. In fact, Tuna and Kirca's solution can be easily verified by plotting both sides of the integro-differential governing equation directly (Fig. 1 in Tuna and Kirca's reply^[24]), which would show that their solutions are valid only for the inner part of the beam.

Wang et al.^[25] obtained exact solutions for the static bending of Euler-Bernoulli beams using Eringen's two-phase local/nonlocal model, where the Fredholm integral equation of the second kind is transformed to its equivalent differential form. The same softening effect was observed.

Eringen's pure nonlocal and two-phase local/nonlocal integral models mentioned above are of strain-driven type^[26]. A singular behavior would appear if one tries to approach the pure nonlocal state from the two-phase local/nolcal state^[27-29]. Romano and Barretta^[26] later proposed the stress-driven nonlocal integral model to overcome the ill-posedness of the strain-driven pure nonlocal integral model. This theory was then extended to a torsional beam problem^[30-31], a thermoelastic beam problem^[32], and a three-dimensional (3D) continua problem^[29]. More interestingly, results from the stress-driven model showed a stiffening effect, while results from the strain-driven model showed a softening effect^{[26, 33]}.

Despite such extensive work on nonlocal elastic beam problems, little work has been done on the case where the axial force cannot be neglected in the beam using Eringen's integral model. Therefore, the present work tries to investigate the indeterminate problem of an Euler-Bernoulli beam with the bending-induced axial force considered. The basic equations for one-dimensional (1D) Euler-Bernoulli beam problem are derived by degenerating from the 3D nonlocal elastic equations. Then, the special cases of a clamped-clamped beam subject to a concentrated force and a uniformly distributed load are studied in detail, since doubly clamped single-wall carbon nanotubes can also be used as actuators^[2]. Numerical examples are given subsequently revealing non-unique solutions caused by improper boundary conditions. Besides, the influence of internal characteristic length and phase parameters on the deflection and axial force is also discussed for different loadings.

2 Degeneration from 3D to 1D Euler-Bernoulli beam problem

The governing equations of nonlocal beam theory can be obtained in several ways, for example, by using the variational principle as indicated in the work of Fernández-Sáez et al.^[21] and Wang et al.^[25]. Although Li et al.^[18] presented a degeneration from two-dimensional (2D) plane-stress problem, Yao^[34] proved earlier that the 2D plane-stress problem may encounter the risk of losing the single-value properties of deflection and strain. Therefore, in order to avoid this trouble, we provide an alternative derivation, that is, degenerating from the 3D problem to 1D Euler-Bernoulli beam problem. It should be noted that this method is inspired by the work of Cowper^[35], where the basic equations of Timoshenko's beam theory were obtained by integrating the equations of 3D elastic theory.

The static 3D equations of nonlocal elasticity are given as follows.

The equilibrium equation is

(1)

where t_ij is the stress tensor of second order, and x_i is the position vector of the particle.

The small deformation equation is

(2)

where ε_ij is the strain tensor of second order, and u_i is the displacement vector.

Eringen's two-phase local/nonlocal constitutive equation for isotropic materials is

(3)

in which

(4)

(5)

where ϖ₁ and ϖ₂ are two phase parameters describing the volume fraction of local and nonlocal effects, respectively, σ_ij is the classic Cauchy stress tensor corresponding to classical elasticity, κ(x-y) is the nonlocal kernel function, and λ_L, μ_L, E, and υ are Lamé constants, Young's modulus, and Poisson's ratio, respectively.

The boundary conditions are

(6)

(7)

where p_i is the surface force, γ_i is the residual surface force induced by the nonlocal effect, and u_i is the prescribed displacement.

For Eqs. (1), (2), (6), and (7), one can refer to Yao^[34] and Huang^[36-38], where the residual nonlocal body force and residual nonlocal angular momentum are proven to be zero. These residuals are assumed to be zero directly without demonstration in early works^[11]. As to the two-phase local/nonlocal constitutive equation (3), one can refer to Eringen^[10-11], Altan^[39], and Polizzotto^[19].

With the above 3D equations, one may obtain the equations for an Euler-Bernoulli beam problem by integrating the 3D equations over its cross section. A straight beam is shown in Fig. 1. Suppose that the x₁-coordinate passes through the geometric centroids of cross sections, the cross section is symmetric about the x₁x₂-plane and the x₁x₃-plane, and all the external loads are applied within the x₁x₂-plane.

Integration of Eq. (1) over the cross section yields

(8)

(9)

From Fig. 1,

(10)

(11)

Then, inserting Eqs. (10) and (11) into Eqs. (8) and (9) and using Stokes' theorem, one may have

(12)

(13)

Equations (12) and (13) are actually the force-equilibrium equations of the beam. To derive the moment-equilibrium equation of the beam, one may need to multiply Eq. (1) by x₂ and integrate over the cross section. Then,

(14)

Similar to the treatment of force-equilibrium equations, Eq. (14) can be written as

(15)

Now define

(16)

Then, equivalent expressions of Eqs. (12), (13) and (15) are

(17)

(18)

(19)

According to Eqs. (17)-(19), it is obvious that F₁(x₁) represents x₁-component of the resultant axial force of the cross section, Q₂(x₁) is the shear force, and M₃(x₁) is the bending moment. Besides, with the existence of p₁, contribution of the axial force to bending moment needs to be taken into consideration. Furthermore, the residual surface force γ_i also influences the axial force, shear force, and bending moment. However, since its representation still remains unknown, its influence is usually supposed to be small and can be neglected when compared with that of external loads in the following analysis.

Now consider the geometric relationship. Denote u₁^*(x₁) and w(x₁) by the x₁-displacement and x₂-displacement (or deflection) of the geometric centroid of cross section. Under the assumption of small deflection, the displacement components expressed in terms of w(x₁) and u₁^*(x₁) are

(20)

Inserting Eq. (20) into Eq. (2) gives

(21)

Based on Eq. (21), the constitutive equation (3) degenerates to

(22)

where A(x₁) is the cross section area of the beam, and κ^*(x₁-y₁) is the nonlocal kernel function for the 1D problem. Equation (22) is consistent with that of Khodabakhshi and Reddy^[20].

Inserting Eq. (22) into Eq. (16), the bending moment M₃(x₁) can be given by

(23)

where I₃ is the moment of inertia with respect to the x₃-axis, and

(24)

The following relationship is used to simplify Eq. (23), that is, due to symmetry of the cross section about the x₂x₃-plane,

(25)

Pay attention that Eqs. (24) and (25) are calculated based on the undeformed state.

Obviously, in the case of classic local beam theory, ϖ₁=1, ϖ₂=0, and the nonlocal kernel function κ(x₁-y₁)=0. Then, Eq. (23) reduces to the classical differential equation of deflection,

(26)

The nonlocal kernel function used here is of bio-exponential form,

(27)

which is the same as that of Erigen^[11]. Substituting Eq. (27) into Eq. (23) yields the same integro-differential governing equation of deflection,

(28)

where a₀ is the internal characteristic length.

3 Semi-analytic solution of clamped-clamped beam subject to concentrated force

As shown in Fig. 2(a), a uniform beam of length 2L fixed at its two ends is subject to a concentrated force 2P at its middle. Due to symmetry of the problem, half of the beam is studied, and the unknown x₁-component X₁ of axial force and bending moment X₃ are labeled in Fig. 2(b).

Taking the axial force into consideration, the shear force Q₂(x₁) and the bending moment M₃(x₁) are given by

(29)

Then,

(30)

Therefore, from Eqs. (30) and (19), the contribution of axial force to bending may be expressed as

(31)

The integro-differential governing equation of deflection is the same as Eq. (28), and the four essential boundary conditions are

(32)

Differentiating Eq. (28) with respect to x₁ twice yields

(33)

Eliminating the integral terms from Eqs. (28) and (33), one obtains

(34)

Differentiating Eq. (34) with respect to x₁ once and introducing

(35)

result in the dimensionless differential equation of rotation angle θ of the cross section,

(36)

where

(37)

The general solution to Eq. (36) is given by

(38)

Then, the deflection of the beam can be derived from Eq. (35),

(39)

where

(40)

Pay attention that the term b₂²-4b₁ is nonnegative, because

(41)

Up till now, there are eight unknown variables (B₁, B₂, B₃, B₄, B₅, w_B, X₁, and X₃) to be determined, which means that 8 boundary conditions are needed. The first 4 obvious boundary conditions are Eq. (32) and can be equivalently written as

(42)

The 5th and 6th conditions may be derived by transferring Eq. (28) into a Fredholm integral equation of the second kind, that is,

(43)

where

(44)

(45)

According to Polyanin and Manzhirov^[40], the function (ξ) should satisfy the following conditions:

(46)

(47)

which may result in

(48)

(49)

Notice that the above 6 boundary conditions (42), (48), and (49) are compulsory, while the following 7th and 8th boundary conditions may be chosen in different ways. For example, according to Eq. (28), the 7th and 8th boundary conditions might be expressed as

(50)

(51)

However, if Eqs. (50) and (51) are both used, numerical calculations in Subsection 6.2 would show that there exist several different solutions satisfying all the 8 boundary conditions and the original integro-differential governing equation (28). Thus, a unique solution cannot be determined. This is quite an interesting phenomenon. Mechanism for causing multi-solutions is indeed very complicated, and much more work needs to be done in the future. Nonlinearity and boundary conditions are two of the many factors that happen to be found in this work.

In this paper, the 7th boundary condition is chosen as Eq. (51), while the 8th boundary equation is found by calculating the extension of the centroid line. Suppose that the displacement u₁ (s₁) of the centroid along the centroid line is a linear function of s₁, which means

(52)

where s₁ is the coordinate along the centroid line, and ε₀ is a constant. According to Eqs. (21) and (22), one may similarly find the strain and stress along the centroid line being

(53)

(54)

where one important property of the kernel function is used to simplify Eq. (54), that is,

(55)

Equation (54) implies the axial force F_N to be constant along the centroid line. Thus, the extension of centroid line can be given by

(56)

where Eq. (42) and the following relation are used,

(57)

Moreover, the extension of centroid line can also be calculated as

(58)

Equating Eqs. (56) and (58), one finally gets the 8th boundary condition,

(59)

Here, we want to make some supplements to the 7th and 8th boundary conditions before we proceed on.

In order to find the solution to the original integro-differential governing equation of deflection (28), we turn to find the solution to its differential equation (34). Therefore, the integro-differential equation (28) can be used as the boundary condition for the differential equation (34). According to Tuna and Kirca^{[22, 24]}, Romano and Berretta^[23], and Wang et al.^[25], special attention needs to be paid to the two ends of the beam. Thus, we choose the 7th boundary condition either as Eq. (50) or as Eq. (51).

As to the 8th boundary condition, it is chosen due to its physical meaning. One might see that the assumption (52) is quite an important one, which results in the constant axial force along the centroid line. Hence, the semi-analytical solution presented in this paper is approximate in this sense.

Now, we have 8 sufficient boundary conditions, which are Eqs. (42), (51), (56), and (59). It is not difficult to see that B₁, B₂, B₃, B₄, B₅, and w_B can be expressed in terms of X₁ and X₃ analytically if one solves Eqs. (42) and (51), since they are liner equations of B₁, B₂, B₃, B₄, B₅, and w_B. For simplicity, the redundant expressions are omitted here. After taking the coefficients into Eqs. (56) and (59), one should arrive at two nonlinear equations about X₁ and X₃. That is why the solution is called semi-analytical. Numerical examples are given in Subsection 6.1 to show that the solution obtained satisfies the original integro-differential governing equation of deflection (28).

For a better numerical calculation when using MAPLE, the above equations (28), (38), (39), (42), (51), (56), and (59) should be transformed into dimensionless ones to avoid extremely large or small values. This can be achieved by introducing

(60)

Then, the equivalent dimensionless expressions are as follows.

For the original integro-differential governing equation,

(61)

where

(62)

For general solutions of θ(ξ) and W(ξ),

(63)

(64)

where

(65)

For the 8 boundary conditions,

(66)

(67)

(68)

(69)

(70)

where

(71)

In Eqs. (61)-(71), the coefficients to be determined are B₁, B₂, ..., B₅, W_B, α₁, and α₃.

4 Semi-analytic solution of clamped-clamped beam subject to uniformly distributed load

If the beam is subject to a uniform distributed load q₀ as shown in Fig. 3, the shear force Q₂(x₁), the bending moment M₃(x₁), and the integro-differential governing equation of deflection are given by

(72)

(73)

(74)

Similar to the process of Section 3, one may obtain the differential form of Eq. (74),

(75)

Differentiating Eq. (75) with respect to x₁ and introducing

(76)

result in

(77)

where

(78)

Besides, the integro-differential equation (75) can be written as

(79)

where

(80)

The general solutions of θ(ξ) and W(ξ) for Eqs. (76) and (75) are given by

(81)

(82)

where

(83)

Similarly, the required first 4 boundary conditions are

(84)

Substituting M₃(x₁) in Eq. (45) with Eq. (73), and taking the updated f(ξ) into Eqs. (46) and (47), one may yield the 5th and 6th boundary conditions,

(85)

(86)

The 7th and 8th boundary conditions are

(87)

(88)

where

(89)

The coefficients H₁, H₂, H₃, H₄, H₅, and W_B can be expressed in terms of α₁ and α₃ by solving Eqs. (84)-(86). Then, the values of α₁ and α₃ can be obtained from nonlinear equations (87) and (88) numerically. Examples are given in Subsection 6.2.

5 Solutions to problems using classical local beam theory 5.1 Case of concentrated force

In the classical local beam theory, that is, without considering the nonlocal size effect, the solution to the problem can be similarly obtained. In this case, ϖ₁=1, ϖ₂=0, and λ₀→ −∞. Then, Eqs. (61) and (36) reduce to

(90)

(91)

where

(92)

The general solutions of θ(ξ) and W(ξ) for Eqs. (90) and (91) are then given by

(93)

(94)

There are 6 unknowns (C₁, C₂, C₃, W_B, α₁, and α₃) to be determined, demanding 6 boundary conditions, which are Eq. (66), Eq. (70), and the following:

(95)

5.2 Case of uniformly distributed load

In this case, ϖ₁=1, ϖ₂=0, and λ₀→ −∞. Then, Eqs. (79) and (77) reduce to

(96)

(97)

where

(98)

The general solutions of θ(ξ) and W(ξ) for Eqs. (96) and (97) are then given by

(99)

(100)

All the 6 unknowns (K₁, K₂, K₃, W_B, α₁, and α₃) can be solved in a similar way, where the required 6 boundary conditions are Eq. (84), Eq. (88), and the following:

(101)

6 Numerical examples and discussion 6.1 Numerical verification of solutions

Consider a nanobeam of length L=100 nm, whose cross section is circular and of diameter d₀=5 nm. Given that Young's modulus E=0.72 TPa^[41] and ϖ₁=1/2. Besides, the concentrated force P =1.0×10^-9 N, and the distributed load q₀=P/L. By using MAPLE, the solutions can be verified numerically by plotting both sides of the integro-differential governing equation (28) or its dimensionless form (61) for concentrated loading (see Fig. 4(a)) and Eq. (79) for distributed loading (see Fig. 5(a)). The corresponding differences between the left side f_left and the right side f_right are given in Figs. 4(b) and 5(b), respectively.

According to Figs. 4 and 5, it is shown that the obtained solutions are of high accuracy, within a tolerance of about 10^-9.

6.2 Illustration of non-unique solutions

It has been mentioned in Section 3 that an improper choice of additional boundary conditions might cause non-unique solutions. Figure 6 illustrates this point, where L=100 nm, d₀= 5 nm, E=0.72 TPa, a₀/L=0.05, ϖ₁=1/2, and q₀=1.0×10^-11 N/nm, and Eqs. (50) and (51) are used as the additional 7th and 8th boundary conditions. Interestingly, nonlocal solution 1 shows a softening effect, while nonlocal solutions 2 and 3 manifest a stiffening effect when compared with the local solution. However, if the 7th boundary condition (50) is replaced by Eq. (88), the unique solution is found to be nonlocal solution 1. Notice that only three nonlocal solutions are displayed in Fig. 6. In fact, there are more solutions.

Khodabakhshi and Reddy^[20] once used the FEM formulation for beam element to solve a clamped-clamped beam with the uniformly distributed load. Whether they considered the axial force or not is uncertain, since the reference deflection they used is actually for the case without the axial force (one can refer to Subsection 3.2 of Ref. [20]). Regardless of this, one can still find that they pointed out one thing that the FEM solutions for the whole beam model and the half beam model due to symmetry had a noticeable difference (one can refer to Figs. 2 and 4 of Ref. [20]). Except for the reasons already given in Ref. [20], we guess that such difference might have something to do with the choice of boundary conditions. Perhaps, they obtained non-unique solutions.

6.3 Influence of internal characteristic length and phase parameters on deflection and axial force of beam

Figure 7 displays the influence of internal characteristic length on the deflection and rotation angle of the cross section for the concentrated loading, while Fig. 8 for distributed loading. Both figures show a consistent softening effect, that is, the deflection and the rotation angle increase as the internal characteristic length increases.

Figure 9 plots the ratio of axial force X₁ to the applied load P or q₀L. For both loading cases, when the load increases, the axial force increases, however, with a decreasing gradient. At a given load, a larger internal characteristic length results in a larger axial force. Besides, for a clamped-clamped beam, the axial force is comparable to or several times as large as the applied load. Therefore, the axial force may not be neglected, since it would also affect the fundamental frequency of nanotubes and nanobeams.

Furthermore, Fig. 10 tells that the maximum deflection w_B decreases as the phase parameter ϖ₁ increases for concentrated and distributed loads. The same trend happens to the axial force, which is not plotted here, since a larger maximum deflection relates to a larger axial force.

7 Conclusions

This paper studies the Euler-Bernoulli beam problem using Eringen's two-phase local/ nonlocal integral constitutive equation. Particularly, a clamped-clamped beam subject to a concentrated force and a uniformly distributed load is solved semi-analytically, where the bending induced axial force is considered. According to the numerical results, we arrive at the following conclusions.

(ⅰ) The axial force is significant and should be taken into consideration in a clamped-clamped beam, and a calculation of shear force and bending moment needs to be based on the deformed state.

(ⅱ) Since for nonlinear problems, an inadequate choice of boundary conditions may give non-unique solutions, boundary conditions should be chosen carefully. By introducing an alternative physically-meaningful boundary condition, a unique solution can be obtained. The method presented here can be applied to other beam problems.

(ⅲ) A consistent softening size effect exists for a clamped-clamped beam using Eringen's two-phase local/nonlocal integral model.

Eringen's two-phase local/nonlocal constitutive equation is a Fredholm equation of the second kind, while Eringen's pure integral constitutive equation belongs to the Fredholm equation of the first kind, which is much more difficult to solve. Solutions in this paper cannot be degenerated to those using Eringen's pure integral model by letting the phase parameter ϖ₁=0, since the two integral equations are different in essence. Much more work needs to be done on Eringen's pure integral model.

The presented semi-analytical results in this paper are aimed to provide useful guidelines for further study of dynamical properties of actuators^[2] or sensors using a clamped-clamped nanobeam structure.