A Nonlocal Model for Carbon Nanotubes under Axial Loads

Abstract

Various beam theories are formulated in literature using the nonlocal differential constitutive relation proposed by Eringen. A new variational framework is derived in the present paper by following a consistent thermodynamic approach based on a nonlocal constitutive law of gradient-type. Contrary to the results obtained by Eringen, the new model exhibits the nonlocality effect also for constant axial load distributions. The treatment can be adopted to get new benchmarks for numerical analyses.

1. Introduction

Carbon nanotubes (CNTs) are a topic of major interest both from theoretical and applicative points of view. This subject is widely investigated in literature to describe small-scale effects [1–4], vibration and buckling [5–13], and nonlocal finite element analysis [14–18]. A comprehensive review on applications of nonlocal elastic models for CNTs is reported in [19] and therein references. Buckling of triple-walled CNTs under temperature fields is dealt with in [20]. An alternative methodology is based on an atomistic-based approach [21] which predicts the positions of atoms in terms of interactive forces and boundary conditions. The standard approach to analyze CNTs under axial loads consists in solving an inhomogeneous second-order ordinary differential equation providing the axial displacement field, see, for example, [22]. The known term of the differential equation is the sum of two contributions. The former describes the local effects linearly depending on the axial load. The latter characterizes the small-scale effects depending linearly on the second derivative along the rod axis of the axial load. This model is thus not able to evaluate small-scale effects due to, constant axial loads per unit length. This approach, commonly adopted in literature, is based on the following nonlocal linearly elastic constitutive law proposed by Eringen [23]: 𝜎−𝑒²_𝑜𝑎²𝜎⁽²⁾=𝐸𝜀, (1) where 𝑒_𝑜 is a material constant, 𝑎 is the internal length, 𝐸 is the Young modulus, 𝜎 is the normal stress, the apex (•)⁽²⁾ is second derivative along the rod axis, and 𝜀 is the axial elongation. Indeed, integrating on the rod cross section domain Ω and imposing that the axial force 𝑁 is equal to the resultant of normal stress field we get the differential equation 𝑁−𝑒²_𝑜𝑎²𝑁⁽²⁾=𝐸𝐴𝑤⁽¹⁾, (2) where 𝜀=𝑤⁽¹⁾, with 𝑤⁽¹⁾ being first derivative along the rod axis of the axial displacement field 𝑤∶[0,𝐿]↦ℛ, where 𝐿 is the rod length and 𝐴 denotes the cross section area. Since the equilibrium prescribes that the first derivative of 𝑁 is opposite to the axial load 𝑝, we infer the well-known differential equation (see, e.g., [7]) as follows: 𝑤⁽²⁾=−𝑝+𝑒²_𝑜𝑎²𝑝⁽²⁾.𝐸𝐴 (3) Note that the nonlocal contribution vanishes for constant loads 𝑝. In the present paper, an alternative nonlocal constitutive behavior is adopted to assess small-scale effects in nanotubes also for constant axial loads. The corresponding axial displacement field is shown to be governed by a fourth-order inhomogeneous differential equation. Boundary conditions are naturally inferred by performing a standard localization procedure of a variational problem formulated by making recourse to thermodynamic restrictions see, for example, [24–26], according to the geometric approach illustrated in [27–30]. As an example, the displacement field of nanotubes under constant axial loads per unit length is evaluated in the appendix. Vibration and buckling effects are not the subject of this paper and will be addressed in a forthcoming paper.

2. Nonlocal Variational Formulation

Let ℬ be the three-dimensional spatial domain of a straight rod subjected to axial loads. An apex (•)^(𝑛) stands for 𝑛th derivative along the rod centroidal 𝑧-axis. Kinematic compatibility between axial elongations 𝜀 and axial displacements 𝑤 is expressed by the differential equation 𝜀=𝑤⁽¹⁾. Denoting by a dot the time-rate, the following noteworthy relations hold true 𝜀̇=𝑤̇⁽¹⁾,𝜀̇⁽¹⁾=𝑤̇⁽²⁾. (4) The differential equation of equilibrium turns out to be 𝑁⁽¹⁾=−𝑝. Boundary equilibrium prescribes that at the end cross sections act axial loads equal to 𝑁(0) for 𝑧=0 and to 𝑁(𝐿) for 𝑧=𝐿. Let us now consider a nonlocal constitutive model of gradient-type defined by assigning the following elastic energy functional per unit volume: 𝜓𝜀,𝜀⁽¹⁾1∶=2𝐸𝜀²+12𝐸𝑐²𝜀²⁽¹⁾, (5) with 𝑐∶=𝑒_𝑜𝑎 being nonlocal parameter. Relation (5) is similar to the elastic energy density proposed in [31] where a homogeneous quadratic functional including also mixed terms is assumed. The elastic energy time rate is hence expressed by the formula 𝜓̇𝜀,𝜀⁽¹⁾=𝜕𝜓̇𝜕𝜀𝜀+𝜕𝜓𝜕𝜀⁽¹⁾𝜀̇⁽¹⁾=𝜎_𝑜𝜀̇+𝜎₁𝜀̇⁽¹⁾, (6) where 𝜎_𝑜∶=𝜕𝜓𝜕𝜀=𝐸𝜀,𝜎₁∶=𝜕𝜓𝜕𝜀⁽¹⁾=𝐸𝑐²𝜀⁽¹⁾, (7) are the static variables conjugating with the kinematic variables 𝜀 and 𝜀⁽¹⁾. The static variable 𝜎₁ is the scalar counterpart of the so-called double stress tensor [31]. By imposing the thermodynamic condition (see, e.g., [32–34]) _ℬ𝜎𝜀𝑑̇𝑉−_ℬ𝜓̇𝑑𝑉=0, (8) where 𝜎 is the normal stress, we infer the relation _ℬ𝜎𝜀𝑑̇𝑉=_ℬ𝜎_𝑜𝜀𝑑̇𝑉+_ℬ𝜎₁𝜀̇⁽¹⁾𝑑𝑉. (9) The relevant differential and boundary equations are thus obtained as shown hereafter. Substituting the expression of the rates 𝜀̇ and 𝜀̇⁽¹⁾ in terms of the axial displacement 𝑤(𝑧) of the cross section at abscissa 𝑧, we get the formulae _ℬ𝜎𝜀𝑑̇𝑉=_ℬ𝜎𝑤̇⁽¹⁾𝑑𝑉=^𝐿_0Ω𝜎𝑑𝐴𝑤̇⁽¹⁾𝑑𝑧=^𝐿₀𝑁𝑤̇⁽¹⁾𝑑𝑧,_ℬ𝜎_𝑜𝜀𝑑̇𝑉=^𝐿₀𝑁_𝑜𝑤̇⁽¹⁾𝑑𝑧,_ℬ𝜎₁𝜀̇⁽¹⁾𝑑𝑉=_ℬ𝜎₁𝑤̇⁽²⁾𝑑𝑉=^𝐿_0Ω𝜎₁𝑑𝐴𝑤̇⁽²⁾=𝑑𝑧_𝐿0𝑁₁𝑤̇⁽²⁾𝑑𝑧, (10) with 𝑁=∫_Ω𝜎𝑑𝐴 axial force (static equivalence condition on the cross sections) and 𝑁_𝑖∶=∫_Ω𝜎_𝑖𝑑𝐴 for 𝑖∈{0,1}. Thermodynamic condition (9) provides the axial contribution ^𝐿₀𝑁𝑤̇⁽¹⁾𝑑𝑧=^𝐿₀𝑁_𝑜𝑤̇⁽¹⁾𝑑𝑧+^𝐿₀𝑁₁𝑤̇⁽²⁾𝑑𝑧. (11)

3. Differential and Boundary Equations of Elastic Equilibrium

Resorting to Green’s formula, a standard localization procedure provides the differential and boundary equations corresponding to the variational conditions inferred in Section 2, as follows. A direct computation gives ^𝐿₀𝑁𝑤̇⁽¹⁾𝑑𝑧=[𝑁𝑤̇]^𝐿₀−^𝐿₀𝑁⁽¹⁾𝑤̇𝑑𝑧,_𝐿0𝑁_𝑜𝑤̇⁽¹⁾𝑑𝑧=𝑁_𝑜𝑤̇^𝐿₀−^𝐿₀𝑁_𝑜⁽¹⁾𝑤̇𝑑𝑧,_𝐿0𝑁₁𝑤̇⁽²⁾𝑑𝑧=𝑁₁𝑤̇^(1)𝐿₀−^𝐿₀𝑁₁⁽¹⁾𝑤̇⁽¹⁾=𝑁𝑑𝑧₁𝑤̇^(1)𝐿₀−𝑁₁⁽¹⁾𝑤̇^𝐿₀+^𝐿₀𝑁₁⁽²⁾𝑤̇𝑑𝑧. (12)

Substituting into the variational condition (11), a suitable localization provides the relevant differential equation 𝑁⁽¹⁾=𝑁_𝑜⁽¹⁾−𝑁₁⁽²⁾ (13) and boundary conditions 𝑁=𝑁_𝑜−𝑁₁⁽¹⁾,dualof𝑤̇,0=𝑁₁,dualof𝑤̇⁽¹⁾. (14)

These conditions can be conveniently expressed in terms of the axial displacement field 𝑤 as follows. A direct evaluation of the scalar functions 𝑁_𝑖∶[0,𝐿]↦ℛ for 𝑖∈{0,1} and of their derivatives gives 𝑁_𝑜=_Ω𝜎_𝑜𝑑𝐴=_Ω𝐸𝜀𝑑𝐴=_Ω𝐸𝑤⁽¹⁾𝑑𝐴=𝐸𝐴𝑤⁽¹⁾,𝑁_𝑜^(𝑗)=𝐸𝐴𝑤^(1+𝑗),𝑁₁=_Ω𝜎₁𝑑𝐴=_Ω𝐸𝑐²𝜀⁽¹⁾=𝑑𝐴_Ω𝐸𝑐²𝑤⁽²⁾𝑑𝐴=𝐸𝐴𝑐²𝑤⁽²⁾,𝑁₁^(𝑗)=𝐸𝐴𝑐²𝑤^(2+𝑗), (15) with 𝑗∈{1,2,…,𝑛}. Accordingly, the boundary and differential conditions of elastic equilibrium (13) and (14) take the form 𝑁⁽¹⁾=𝑁_𝑜⁽¹⁾−𝑁₁⁽²⁾=𝐸𝐴𝑤⁽²⁾−𝐸𝐴𝑐²𝑤⁽⁴⁾,𝑁=𝑁_𝑜−𝑁₁⁽¹⁾=𝐸𝐴𝑤⁽¹⁾−𝐸𝐴𝑐²𝑤⁽³⁾,dualof𝑤̇,0=𝑁₁=𝐸𝐴𝑐²𝑤⁽²⁾,dualof𝑤̇⁽¹⁾. (16)

4. Example

Let us consider a straight rod subject to a constant axial load 𝑝 as depicted in Figure 1. End cross sections 𝐀 and 𝐁 are assumed to be hinged and simply supported, respectively. As illustrated in Section 3, the computation of the rod axial displacement field 𝑤 involves the following cross section geometric and elastic properties: area 𝐴, Young modulus 𝐸, and nonlocal parameter 𝑐. By setting 𝛼∶=𝐸𝐴𝑐² and 𝛽∶=𝐸𝐴, the differential equation of elastic equilibrium is as follows: 𝛼𝑤⁽⁴⁾−𝛽𝑤⁽²⁾=𝑝. (17) The general integral takes thus the form (see the appendix) 𝑤(𝑧)=𝑤_𝐻(𝑧)+𝑤(𝑧), (18) with 𝑤_𝐻(𝑧)=𝑐₁+𝑐₂𝑧+𝑐₃exp𝛽𝛼𝑧+𝑐₄exp−𝛽𝛼𝑧,𝑤(𝑧)=𝑐₅𝑧². (19) The evaluation of the constants is carried out by imposing the following boundary conditions (see also Section 3): 𝑤(0)𝑁(𝐿)=0,=𝐸𝐴𝑤⁽¹⁾(𝐿)−𝐸𝐴𝑐²𝑤⁽³⁾(𝐿)𝑁=0,₁(0)=𝐸𝐴𝑐²𝑤⁽²⁾(0)𝑁=0,₁(𝐿)=𝐸𝐴𝑐²𝑤⁽²⁾(𝐿)⇕𝑤(0)𝑤=0,=0,₍₁₎(𝐿)−𝑐²𝑤⁽³⁾(𝐿)𝑤=0,₍₂₎(0)𝑤=0,₍₂₎(𝐿)=0. (20) Resorting to the expressions of the derivatives 𝑤_𝐻^(𝑗) and 𝑤^(𝑗) for 𝑗∈{1,2,3,4}, 𝑤_𝐻⁽¹⁾(𝑧)=𝑐₂+𝑐₃𝑐1exp𝑐𝑧−𝑐₄𝑐−1exp𝑐𝑧,𝑤_𝐻⁽²⁾(𝑧)=𝑐₃𝑐²1exp𝑐𝑧+𝑐₄𝑐²−1exp𝑐𝑧,𝑤_𝐻⁽³⁾(𝑧)=𝑐₃𝑐³1exp𝑐𝑧−𝑐₄𝑐³−1exp𝑐𝑧,𝑤_𝐻⁽⁴⁾(𝑧)=𝑐₃𝑐⁴1exp𝑐𝑧+𝑐₄𝑐⁴−1exp𝑐𝑧,𝑤⁽¹⁾(𝑧)=2𝑐₅𝑧,𝑤⁽²⁾(𝑧)=2𝑐₅,𝑤⁽³⁾(𝑧)=0,𝑤⁽⁴⁾(𝑧)=0, (21) and having 𝛽/𝛼=1/𝑐, a direct computation provides the algebraic system 𝑐₁+𝑐₃+𝑐₄𝑐=0,₂+2𝐿𝑐₅1=0,𝑐²𝑐₃+1𝑐²𝑐₄+2𝑐₅1=0,𝑐²1exp𝑐𝐿𝑐₃+1𝑐²−1exp𝑐𝐿𝑐₄+2𝑐₅=0. (22) A further condition can be obtained by imposing that the scalar field 𝑤(𝑧)=𝑐₅𝑧² (23) is a particular solution of the differential equation (17), whence it follows that 𝑐₅=−𝑝/2𝐸𝐴. The remaining constants are given by the formulae 𝑐₁=−𝑐₃+𝑐₄,𝑐₂=−2𝐿𝑐₅,𝑐₃=2𝑐²1−exp(−(1/𝑐)𝐿)(−()𝐿))𝐿)𝑐exp1/𝑐−exp((1/𝑐⁵,𝑐₄=2𝑐²1−exp((1/𝑐)𝐿))𝐿)(−()𝐿)𝑐exp((1/𝑐−exp1/𝑐⁵, (24) having 𝑐₅=−𝑝/2𝐸𝐴. A plot of the rod axial displacement field 𝑤 for different values of the nonlocal parameter 𝑐 is provided in Figure 2. It is apparent that the rod becomes stiffer if the nonlocal parameter increases. The evaluated axial displacement at the free end of the rod 𝐁 provides the same value independently of the nonlocal parameter. Such a value coincides with the displacement of the point 𝐁 if a local model is considered. Moreover, the limit of the axial displacement field for 𝑐 tending to plus infinity can be evaluated to get the lower bound 𝑤^low(𝑧)∶=lim_𝑐→+∞𝑤(𝑧,𝑐)=0.208333𝑧. (25) Hence, large values of the nonlocal parameter provide a displacement field which tends to a linear one, see Figure 2, for 𝑐=25. Further, the limit value of the axial displacement for 𝑧=𝐿 and 𝑐→+∞, obtained by (25), yields 𝑤^low(𝐿)=2.08333 nm which coincides with the axial displacement at 𝐁 for any value of the nonlocal parameter 𝑐, see Figure 3 and Table 1.

Table 1 

Axial displacements 𝑤(𝐿/2) and 𝑤(𝐿) versus the nonlocal parameter 𝑐.

Figure 1 

Rod under constant axial load.

Figure 2 

Axial displacement field 𝑤 for 𝑐=0 (local solution—black line), 𝑐=1 nm (blue line), 𝑐=2 nm (red line), 𝑐=3 nm (green line), 𝑐=4 nm (orange line), 𝑐=5 nm (yellow line), and 𝑐=25 nm (brown line). 𝐸=300 GPa; 𝐿=10 nm; 𝐴=80⋅10−2 nm2; and 𝑝=10−8 N/nm.

Figure 3 

Axial displacement in terms of the nonlocal parameter 𝑐 at the abscissa 𝑧=𝐿/2 (blue line) and 𝑧=𝐿 (red line).

The upper bound of the axial displacement is provided by the local solution (i.e., 𝑐=0) 𝑤^upp(𝑧)()=𝑝(𝑧)𝑧∶=𝑤𝑧,0)(2𝐿−𝑧.2𝐸𝐴 (26) The axial displacement evaluated for 𝑧=𝐿 by (26) yields the value 𝑤^upp(𝐿)=25/12 nm which coincides with the axial displacement at 𝐁 for any value of 𝑐, see Figure 3 and Table 1. For the considered model, the upper and lower bounds of the axial displacement field are given by (25) and (26). The axial displacement 𝑤(𝐿/2) at the middle point of the rod and the maximum axial displacement 𝑤(𝐿) as functions of the nonlocal parameter 𝑐 are depicted in Figure 3. The corresponding numerical values of 𝑤(𝐿/2) and 𝑤(𝐿) are listed in Table 1.

It is worth noting that equilibrium prescribes that axial force 𝑁 must be a linear function, confirmed by the blue diagram in Figure 4 obtained as difference between the local contribution 𝑁_𝑜 (dashed line) and the nonlocal one 𝑁₁⁽¹⁾ (continuous thin line), according to (14)₁ for any value of 𝑐.

Figure 4 

Axial force 𝑁=𝑁𝑜−𝑁1(1) (blue line), 𝑁𝑜 (dashed line), 𝑁1(1) (continuous thin line); 𝑐=0 (local solution-black line), 𝑐=1 nm (blue line), 𝑐=2 nm (red line), 𝑐=3 nm (green line), 𝑐=4 nm (orange line), 𝑐=5 nm (yellow line), and 𝑐=25 nm (brown line).

5. Conclusions

The outcomes of the present paper may be summarized as follows. (i) Linearly elastic carbon nanotubes under axial loads have been investigated by a nonlocal variational approach based on thermodynamic restrictions. The treatment provides an effective tool to evaluate small-scale effects in nanotubes subject also to constant axial loads, a goal not achievable by the Eringen model commonly adopted in literature as motivated in Section 1. (ii)Relevant boundary and differential conditions of elastic equilibrium have been inferred by a standard localization procedure. Such a procedure provides, in a consistent way, the relevant class of boundary conditions for the nonlocal model. (iii)The present approach yields a firm thermodynamic procedure to derive different nonlocal models for CNTs by suitable specializations of the elastic energy. (iv)Exact solutions of carbon nanotubes subject to a constant axial load have been obtained. An advantage of the proposed procedure consists in providing an effective tool to be used as a benchmark for numerical analyses. Finally, a range to which any nonlocal solution must belong is analytically evaluated.

Appendix

The procedure to solve the ordinary differential equation 𝛼𝑤⁽⁴⁾−𝛽𝑤⁽²⁾=𝑓, (A.1) with 𝛼, 𝛽>0 being constant coefficients and 𝑓∶𝐼⊆ℛ↦ℛ being a continuous function, is summarized as follows. Let us consider the homogeneous differential equation 𝛼𝑤⁽⁴⁾−𝛽𝑤⁽²⁾=0 (A.2) and the relevant characteristic (algebraic) equation 𝛼𝜆⁴−𝛽𝜆²=0. The roots of the polynomial 𝛼𝜆⁴−𝛽𝜆² are 𝜆₁=0 with multiplicity 2, 𝜆₂=𝛽/𝛼 with multiplicity 1 and 𝜆₃=−𝛽/𝛼 with multiplicity 1. The general integral of (A.2) is thus expressed by the formula 𝑤_𝐻(𝑧)=𝑐₁+𝑐₂𝑧+𝑐₃exp𝛽𝛼𝑧+𝑐₄exp−𝛽𝛼𝑧, (A.3) with exp denoting exponential function and 𝑐_𝑖∈ℛ for 𝑖={1,…,4}. The general integral of (A.1) is writen therefore as 𝑤(𝑧)=𝑤_𝐻(𝑧)+𝑤(𝑧), (A.4) where 𝑤 is a particular solution of (A.1). It is worth noting that, for 𝑓 defined by a polynomial 𝑝_𝑚 of degree 𝑚≥0, the solution 𝑣 can be looked for by setting 𝑤(𝑧)=𝑧²𝐴₀+𝐴₁𝑧+⋯+𝐴_𝑚𝑧^𝑚, (A.5) with 𝐴_𝑖∈ℛ for 𝑖={1,…,𝑚}.

Acknowledgments

The authors were supported by the “Polo delle Scienze e delle Tecnologie,” University of Naples Federico II, through the research project FARO. Useful hints and precious comments by anonymous reviewers are also gratefully acknowledged.