- Email: [email protected]
On functionally graded Timoshenko nonisothermal nanobeams
Accepted Manuscript On Functionally Graded Timoshenko Nonisothermal Nanobeams Marko Čanađija, Raffaele Barretta, Francesco Marotti de Sciarra PII: DOI: Reference:
S0263-8223(15)00873-9 http://dx.doi.org/10.1016/j.compstruct.2015.09.030 COST 6875
To appear in:
Composite Structures
Please cite this article as: Čanađija, M., Barretta, R., de Sciarra, F.M., On Functionally Graded Timoshenko Nonisothermal Nanobeams, Composite Structures (2015), doi: http://dx.doi.org/10.1016/j.compstruct.2015.09.030
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
On Functionally Graded Timoshenko Nonisothermal Nanobeams Marko Čanađijaa,∗, Raffaele Barrettab , Francesco Marotti de Sciarrab a
Faculty of Engineering,Department of Engineering Mechanics,University of Rijeka,Vukovarska 58,51000 Rijeka,Croatia b Department of Structures for Engineering and Architecture,University of Naples Federico II,Via Claudio 21,80121 Naples,Italy
Abstract The paper deals with the mechanics of nanobeams within nonisothermal environments. Mechanical behaviour is assumed to be time-independent and suitable for the statical beam deformation. The nanobeam mechanics is based on the Timoshenko kinematics, augmented by nonlocality effects as advocated by Eringen. Following a thermodynamic approach, a sound framework oriented towards anisotropic materials is developed. Such nanobeam model is especially suited for functionally graded materials. The proposed procedure is tested on two examples. The first example investigates mechanical behaviour in the more general way for the different temperature distributions, while the second one is oriented toward specific carbon nanotube Poly(methyl methacrylate) nanocomposite. Keywords: Timoshenko nanobeams, nonlocal thermoelasticity, carbon nanotubes 1. Introduction A great interest of researchers is nowadays directed toward carbon nanotubes (CNT) reinforced polymers. This is of course related to the superior CNT mechanical properties, but additional benefits can be further achieved by the functionally grading of the composite [1]. The theoretical research of the subject can be divided into three main directions - static and buckling, ∗
Corresponding author. Tel.: +385-51-651-496; Fax.: +385-51-651-490 Email address: [email protected] (Marko Čanađija )
Preprint submitted to Composite Structures
September 23, 2015
ACCEPTED MANUSCRIPT
dynamic and vibrations and nonlinear analysis. Remarkably, only a tiny fraction of published papers addresses nonisothermal behaviour. Within the isothermal context, previous non-local functionally graded (FG) beam models were developed to address vibration issues [2, 3, 4, 5, 6, 7]. Interests in the vibrational response of a nanobeam lies in the frequent application in nanosenors. In particular, the resonant frequency is exploited in mass sensing at the nanoscale. Another class of sensors are based on the classical beam statics. Hence, bending and buckling of a non-local FG nanobeam was investigated in [8, 9] for Timoshenko formulation and Bernoulli-Euler nanobeam in [10, 11]. Torsion was investigated in [12]. Extension to viscoelastic behaviour was recently proposed in [13]. Beside nanosensorics, application of FG nanocomposites is very wide; it seems that tissue engineering can significantly benefit from such materials, resembling composition and structure of the bone tissue much better than the conventional materials, [14]. FG nanotubes can be used for fluid transport in the nano electromechanical systems [15]. Different grading distributions can be conveniently obtained by the centrifugation techniques [16], resulting in polymer composites with superior wear resistance. Another frequent manufacturing process of FG nanocomposites uses assembling of different layers of materials. For example, aluminium - CNT composite is produced from layers containing different fractions of CNT [17]. In this particular case, Vickers hardness of each particular layer was measured with results varying by a factor of 7. However, when influence of temperature on the nanomechanics is considered, there are only a few papers available. Bernoulli-Euler kinematics of nonisothermal nanobeams for isotropic materials is analyzed in [18]. Previous FG material thermomechanical applications at nanoscale involve bifurcation analysis of functionally graded (FG) shells [19] and nonlinear vibrations of FG shells [20]. Thermal buckling of curved nanowires, yet with a homogeneous temperature field was investigated in [21]. Nonlocality, by means of surface effects were accounted for. It was found that thermal effects significantly affect thermal load by reducing buckling load. Wave propagation in piezoelectric nanobeams, influenced by surface and thermal effects were studied in [22]. The authors pointed out the profound influence of the scale coefficient and temperature on the frequency of aforementioned nanobeams. An extension to temperature dependent properties is provided in [23, 24], with a special emphasis of the transient heating in the form of sinusoidal pulse and ramp-type heating on the vibrations of the nanobeam. It was found that the pulse width significantly affects mechanical response of the 2
ACCEPTED MANUSCRIPT
nanobeam. Other contributions to the vibration analysis of FG nanobeams can be found in [25]. Effect of electrostatic force and thermal shock on capacitive FG nanocantilevers was studied in [26]. Finally, it should be pointed out that thermal properties of FG materials are highly affected by the gradient pattern, resulting in a variety of temperature distributions, [27, 28], thus supporting the need for the reliable non-local thermomechanical model. A standard procedure to predict a material property of a nanocomposite, is based on the classical rule of mixtures. For some property P , there is: P F G = P a V a + P b Vb ,
where Va + Vb = 1.
(1)
Above, indexes a, b denote components a and b, respectively, while V represent the volume fraction. The volume fraction of the component is then specified to vary by some functional dependence, for example linear [5, 29], power law [30, 6, 3], Flory-Rehner equation [31]. Sometimes, CNT properties are scaled in order to better match experimentally noticed behaviour, [5, 19, 29]. The scaling parameter accounts for various defects in the reinforcing material or non-perfect load transfer at the CNT-matrix interface. The same approach is pursued in the present paper as well, although the paper does not set any restrictions on the particular form of functional dependence of a property. The novel aspects of this paper can be stated as follows. Starting from the strict continuum mechanics framework, augmented by the Eringen nonlocal theory [32], nonlocal thermoelastic Timoshenko beam governing equations for the bending moment, axial and shear forces are developed. However, to accommodate different nonlocal effects due to normal and shear stresses as usually assumed in Timoshenko nanobeams [33], an extension of the original formulation proposed by Eringen to the anisotropic form is now proposed. In contrast to other papers in the field, solution is offered for the case of arbitrary temperature field; usual approach assumes a homogeneous temperature field in the whole beam. Such nonisothermal and nonlocal Timoshenko formulation suitable for FG nanocomposite materials with arbitrary temperature fields was not previously addressed in the literature. Additionally, usual approach relies on the formulation of the set of ordinary differential equation (ODE) for the axial and transverse displacements as well as for the rotation angle. Unfortunately, such approach in the problems at hand leads to highly coupled set of ODE and represent an unnecessary complication. We show that much simpler way to solve the problem is determining the 3
ACCEPTED MANUSCRIPT
forces and moments first and subsequent solving for the unknown displacements and remaining integration constants. Such a simplification is most welcomed since benchmarks introduced by means of examples shows that resulting displacement fields provide rather lengthy form even for the simplest cases. The volume of the paper is divided into three main sections. In the second section, a rigorous nonlocal continuum framework is set-up. Usage of such results is performed in the third section to obtain Timoshenko nanobeam governing equations. Example section provides new benchmarks and illustrate main results of the paper. Finally, the conclusion section completes the paper and points out the most significant results. 2. Nonlocal continuum mechanics 2.1. Nonlocal stress tensor In order to introduce small-size effects, we pursue the approach introduced by Eringen [34]. Let us consider the current configuration of a body B ⊂ R at time t. It is assumed that in each point x of the body, the balance of momentum is valid: ¯ + f = ρu, ¨ divσ (2) ¯ is the nonlocal stress tensor at the point x, ρ(x, t) is the mass where σ density, u(x, t) is the displacement vector and f (x, t) are volume forces. According to [32], the nonlocal stress tensor has the form: ¯ σ(x) =
Z
α (|x0 − x| , τNL ) σ(x0 )dV (x0 ).
(3)
B
Above, α (|x0 − x| , τNL ) is the nonlocal modulus and σ(x0 ) is the classical stress tensor at point x0 . Usually, τNL = c/l, the nonlocal parameter c = e0 a, where e0 is a materially dependent constant, a is the internal characteristic length like granular distance or lattice parameter and l is the external characteristic length like crack length or similar. In that way, the stress state at x depends on stress states in all other points x0 ∈ B. Regarding the nonlocal modulus α, it can be shown that α is a Green’s function of a linear differential operator L, i.e. [32]: Lα (|x0 − x| , τNL ) = δ (|x0 − x|) .
4
(4)
ACCEPTED MANUSCRIPT
The above form includes the Dirac delta function δ ensuring that the classical (local) elasticity is encompassed in the theory for the choice c = 0, see [32] for ¯ = σ. Moreover, more details. Application of this result to Eq. (3) gives Lσ in the case when L is a differential operator with constant coefficients, body forces are not significant and inertial loads can be neglected, the classical equation of equilibrium is conveniently obtained, divσ = 0. At this point some constitutive law defining the nonlocal modulus should be introduced. Again, based on the short review of such models in [32], for this paper a suitable choice is: q
˙ 0 − x) (x0 − x)(x 1 K ( ), α (|x0 − x| , c) = 0 2πc2 c
(5)
where K0 is the modified Bessel function. For such a choice, the linear differential operator L can be shown to be [32]: L = 1 − c2 ∇2 ,
(6)
¯ and the classical stress thus providing the relation between the nonlocal σ tensor σ: ¯ = σ. (1 − c2 ∇2 )σ (7) Therefore, the nonlocal stress tensor is influenced by the scalar nonlocal parameter c that leads to the isotropic nonlocal behaviour. At this point we generalize this statement, allowing anisotropic nonlocal behaviour through the new specific tensorial form: h
i
¯ = σ, 1 − c1 ∇2 σ
(8)
¯ = σ + c1 σ1 . σ
(9)
or For the specific choice c1 = c2 I, where I represents Kronecker’s tensor, the classical Eringen form, Eq. (7), is recovered. The anisotropic form is particularly suitable for the nonlocal analysis of Timoshenko nanobeams as it will be shown later. 2.2. Balance equations Now, with the previously introduced nonlocal stress tensor, it is necessary to consider the first and the second law of thermodynamics. These laws can 5
ACCEPTED MANUSCRIPT
be written as:
¯ : d + ρr + P ρe˙ + divq = σ q − div , ρη˙ ≥ ρr θ θ
(10)
where e(x, t) the internal energy, d is the rate of deformation tensor, η is the specific entropy per unit mass, θ is temperature, r(x, t) is the heat supply per unit mass. The heat flux vector across a surface with the normal n is denoted as q. The variable P represents the rate at which the energy is supplied to a point from its neighbourhood. Note that balance of energy at the global R level implies that this rate must vanish, i.e. PdV = 0. Consequently, the Ω
balance of energy can be written in the global form as: Z
(ρe˙ + divq) dV =
Ω
Z
¯ : d + ρr) dV, (σ
(11)
Ω
while the second law of thermodynamics can be rephrased as [33]: R
¯ : ε˙ − σ : ε˙ − σ1 ∗ ∇(1) ε˙ dV = 0 σ
Ω
R Ω
(12)
− qθ grad θ dV ≥ 0.
2.3. Constitutive model The material at hand is assumed to be of the thermoelastic nature, with the nonlocal influence on the stress state in the point. Therefore, the internal energy takes the form e = e(ε, ∇ε, η), where ∇ε reflects nonlocality. At this point we do not make any assumption on the functional dependence of a material property M , so mechanical properties are considered to be given in a general form as functions of three spatial coordinates x = x(x, y, z), where x ∈ B. Regarding thermal behaviour, all properties are considered to be temperature dependent, so M = M (x, θ) = M (x, y, z, θ). For such a case, the Helmholtz free energy function ψ = e − ηθ takes the form: ¯ i ∗ ∇ε − β : ε(θ − θ0 )− ρψ(ε, ∇ε, θ) = 21 ε : C : ε + 12 ∇ε ∗ C 1 − 2θ cε (θ − θ0 )2 ,
(13)
where quantities in these expressions are: C=
∂ 2ψ , ∂ε2
¯i = C
∂ 2ψ , ∂(∇i ε)2
−β =
6
∂ 2ψ , ∂ε∂θ
1 ∂ 2ψ − cε = . θ ∂θ2
(14)
ACCEPTED MANUSCRIPT
Above, cε is the specific heat at constant strain, β is the second order tensor describing anisotropic volumetric expansion due to thermal effects, C is the ¯ i is a higher order fourth order tensor describing local elastic behaviour and C tensor describing elastic effects due to strain gradients. In the case of isotropic volumetric heat expansion, volumetric expansion is reduced to: β = 3καI,
(15)
where κ is the bulk modulus and α is the coefficient of linear thermal expansion (CTE). The anisotropic tensor describing nonlocal parameters introduced in Eq. (8) is assumed to have the specific form: c2 b2 b2 2 c1 = b c2 b2 , b2 b2 c2
(16)
thus assuming the different nonlocal parameters for normal (c) and shear (b) stresses. With the above definition, the following constitutive relationship will be used in this work: ¯ i = c1 · C = c1im Cmjkl . C
(17)
Applied to the classical elasticity, this gives; h
i
¯ i = c2 Iim + b2 (1i 1m − Iim ) [λImj Ikl + µ(Imk Ilj + Iml Ijl )] , C h
(18)
iT
where λ, µ are Lame constants and 1T = 1 1 1 . 2.4. Timoshenko beam kinematics In this section, previously developed continuum mechanics framework is specialized to the nanobeams kinematics operating in the nonisothermal environment. The in-plane bending is assumed to take place in the x − y plane, where x denotes the beam’s longitudinal axis. In order to avoid unnecessary distractions, both mechanical and thermal loading are taken to be independent of the z coordinate. Therefore, in an arbitrary beam point defined by coordinates (x, y, z), the displacement field (ux , uy , uz ) is governed
7
ACCEPTED MANUSCRIPT
by the following equations: ux (x, y) = u0 (x) − ϕ(x)y,
uy (x) = v(x),
uz = 0,
(19)
where u0 (x) represents the axial displacement caused by the axial force and the thermal dilatation originating from the homogeneous part of the temperature field, v is the transverse displacement and the rotation of the beam cross-section is ϕ. Consequently, components of the strain tensor follow in the standard manner: (1)
εx (x, y) = ε = ∂x ux (x, y) = u0 (x) − ϕ(1) (x)y, γxy (x, y) = γ = ∂x uy (x) + ∂y ux (x, y) = v (1) (x) − ϕ(x),
(20)
where exponent (•) implies derivative of the respective order with respect to the longitudinal coordinate. Eq. (20)2 implies that the rotation of the crosssection ϕ(x) is not equal to the rotation of the beam longitudinal axis v (1) (x) but it differs from it for γxy (x, y). Also, we follow the standard practice in the thermomechanics of beams and additively decompose εx into the mechanical εM and the thermal part εT . With the defined Timoshenko nanobeam kinematics framework , the internal energy can be now rephrased as e = e(ε, ε(1) , γ, γ (1) , η). The Helmoltz free energy is then: ψ(ε, ε(1) , γ, γ (1) , θ) = e − ηθ. (21) With the above Legrende transformation and the balance of energy the second law of thermodynamics, Eq. (10)2 takes the form: q σ ¯ : ε˙ + τ¯ : γ˙ − ρψ˙ − ρη θ˙ − gradθ + P ≥ 0, θ
(22)
˙ ∇ε, γ, ∇γ, θ): and introducing the time rate of the Helmholtz energy ψ(ε, ˙ (¯ σ − ρ∂ε ψ) : ε˙ + (¯ τ − ρ∂γ ψ) : γ˙ − ρ (η + ∂θ ψ) θ− q (1) (1) −ρ∂ε(1) ψ ε˙ − ρ∂γ (1) ψ γ˙ − θ gradθ + P ≥ 0.
(23)
The standard arguments [35] along with Eq. (16) are now extended to account for the nonlocal behaviour and lead to the conjugacy between stress and
8
ACCEPTED MANUSCRIPT
strain: σ = ρ∂ε ψ = σM + σT c2 σ1 = ρ∂ε(1) ψ, τ = ρ∂γ ψ, b2 τ1 = ρ∂γ (1) ψ,
(24)
so the second law of thermodynamics, see Eq. (12)1 , is now: R Ω R Ω
σ ¯ ε˙ + τ¯γ˙ − σ ε˙ − τ γ˙ − c2 σ1 ε˙(1) − b2 τ1 γ˙ (1) dV = 0, (25)
− qθ grad θ dV ≥ 0,
where dV = dx da and da = dy dz. At this point, the constitutive behaviour reported in Sec. (2.3) is specialized for the presented Timoshenko beam kinematics. In the fashion of the classical rule of mixture, Eq. (1), we introduce the following form of the free energy function at x ∈ B: ψ(x, ε, ε(1) , θ) = ψa (ε, ε(1) , θ)Va (x) + ψb (ε, ε(1) , θ)Vb (x),
(26)
where ψa , ψb are Helmholtz free energy functions describing two components. Volume fractions of each component Va , Vb are position dependent. The extension to more than two components can be carried out in straightforward manner. Restricting our elaborations to in-plane bending, i.e. dropping dependence on the coordinate z, the following form of the Helmholtz free energy for each component can be used: ρψ(ε, ε(1) , θ) = 21 E(θ)ε2 + 12 c2 E(θ)(ε(1) )2 + 1 χG(θ)γ 2 + 21 b2 χG(θ)(γ (1) )2 − 2 1 cε (θ)(θ − θ0 )2 , −E(θ)α(θ)ε(θ − θ0 ) − 2θ
(27)
where E is the Young’s modulus, G is the shear modulus and χ is the shear correction factor. Alternatively, in the case when the rule of mixture is not an appropriate selection, one can use the more general form: ρψ(x, ε, ε(1) , θ) = 21 E(x, y, θ)ε2 + 12 c2 E(x, y, θ)(ε(1) )2 + 1 χG(x, y, θ)γ 2 + 21 b2 χG(x, y, θ)(γ (1) )2 − 2 1 −E(x, y, θ)α(x, y, θ)ε(θ − θ0 ) − 2θ cε (x, y, θ)(θ − θ0 )2 ,
(28)
instead of Eq. (27). The extension to out-of-plane bending leads to the somewhat more complicated form, but it is straightforward and is left as an exercise to the reader. 9
ACCEPTED MANUSCRIPT
Also, to avoid unnecessary complication in the experimental determination of non-local parameters c and b, these parameters will be treated as temperature independent. 3. Internal forces in Timoshenko beams The starting point for the analysis of internal forces in a Timoshenko beam is Eq. (25)1 . By the application of Timoshenko beam kinematics, Eq. (20), Eq. (25)1 can be transformed into: R
(1)
(1)
(¯ σ (u˙ 0 − ϕ˙ (1) y) − (σM + σT )(u˙ 0 − ϕ˙ (1) y)−
Ω
(29)
(2)
˙ −c2 σ1 (u˙ 0 − ϕ˙ (2) y) + τ¯(v˙ (1) − ϕ)− (1) 2 (2) (1) −τ (v˙ − ϕ) ˙ − b τ1 (v˙ − ϕ˙ ))dv = 0.
In the case of the rule of mixtures, Eq. (1,26,27), the stresses Eq. (24) take the following form: σM = (Va Ea + Vb Eb )ε = EFG ε σT = −(Va Ea αa + Vb Eb αb )(θ − θ0 ) = −(Eα)FG (θ − θ0 ) c2 σ1 = c2 (Va Ea + Vb Eb )ε(1) = c2 EFG ε(1) b2 τ = χ(Va Ga + Vb Gb )γ = GFG γ b2 τ1 = χb2 (Va Ga + Vb Gb )γ = χb2 GFG γ (1) ,
(30)
where effective moduli (•)FG were introduced. Alternatively, if the material properties are defined by Eqs. (27) instead of the rule of mixtures, the appropriate functional dependence must be chosen for (•)FG . Bending moments are introduced as: ¯ = −R σ M ¯ yda, A
(1)
MM = − σM yda = − EFG εyda = ϕ(1) EFG y 2 da − u0 R
R
R
A
A
A R (2)
M1 = − σ1 yda = − EFG ε(1) yda = ϕ R
R
AR
RA
A
A
A
EFG y 2 da −
R
EFG yda,
A (2) R u0 A
EFG yda,
MT = − σT yda = (Eα)FG (θ − θ0 )yda, (31)
10
ACCEPTED MANUSCRIPT
while the shear and the axial forces are: R T¯ = τ¯da,
T =
A R
τ da =
AR
T1 =
R
χGFG γda = v (1) χ GFG da − ϕχ GFG da, R
AR
A
A
A
χGFG γ (1) da = v (2) χ GFG da − ϕ(1) χ GFG da,
τ1 da =
R P¯ = σ ¯ da,
R
A
R
R
A
A
(32)
A
PM =
R
σM da =
A
P1 = PT =
R AR A
R
EFG εda = −ϕ
(1)
A
σ1 da = EFG ε(1) da = −ϕ R
R
EFG yda +
A R (2)
A
(1) u0
EFG yda +
A
R
EFG da,
A (2) R u0 A
EFG da,
σT da = − (Eα)FG (θ − θ0 )da. R
A
Note that due to the fact that the material is functionally graded, the above integrals involve position dependent moduli and CTE. Therefore, these general integrals cannot be simplified to the usual form involving the area and the first and the second moment of inertia. However, to obtain equations in the more familiar form, the following notation for stiffnesses in a cross-section are introduced: kEI =
R
EFG y 2 da kES =
R
EFG yda,
R
AR
A
A
A
kGA = χ GFG da kEA =
EFG da.
(33)
These values are not constants but rather functions of the longitudinal coordinate x. With the above notation, Eqs. (31,32) become now: (1)
MM = ϕ(1) kEI − u0 kES , (2) M1 = ϕ(2) kEI − u0 kES , T = v (1) kGA − ϕkGA = (v (1) − ϕ)kGA , T1 = v (2) kGA − ϕ(1) kGA = (v (2) − ϕ(1) )kGA , (1) PM = −ϕ(1) kES + u0 kEA , (2) P1 = −ϕ(2) kES + u0 kEA .
11
(34)
ACCEPTED MANUSCRIPT
Eq. (29) can be now recast into the form involving bending moments, shear and axial forces: RL
(2)
(1)
[(P¯ − PM − PT )u˙ 0 − c2 P1 u˙ 0 + 0 ¯ − MM − MT + b2 T1 )ϕ˙ (1) − c2 M1 ϕ˙ (2) + +(M (T¯ − T )v˙ (1) − b2 T1 v˙ (2) + (−T¯ + T )ϕ]dx ˙ = 0,
(35)
where L represents the beam length. Application of the integration by parts, for example for the transverse displacement v˙ leads to: RL 0
RL
(•)v˙ (1) dx = [(•)v] ˙ L0 − (•)(1) vdx, ˙ 0
RL
h
(2)
(•)v˙ dx = (•)v˙
(1)
iL
h
− (•) v˙
0
0
(1)
iL 0
RL
(36) (2)
+ (•) vdx. ˙ 0
The same procedure should be applied to the variables u˙ 0 and ϕ as well. In such manner, Eq. (35) takes the following form: h
(P¯ − PM − PT )u˙ 0 h
− c
2
i (1) L P1 u˙ 0 0
h
+ c
iL
2
0
RL
(1)
0 iL (1) P1 u˙ 0 0
RL 0
¯ − MM − MT + b2 T1 )ϕ˙ + (M 2
− c M1 ϕ˙
(1)
h
iL 0
+ (T¯ − T )v˙ h
− b2 T1 v˙ (1)
iL
iL 0
0
h
+ c
2
(2)
− c2 P1 u˙ 0 dx+
h
h
(1)
− (P¯ (1) − PM − PT )u˙ 0 dx−
iL
(1) M1 ϕ˙ 0
iL 0
L
R (1) (1) (1) ¯ (1) − MM ˙ − MT + b2 T1 )ϕdx− − (M
RL
0
(2)
2
− c M1 ϕdx+ ˙ 0
RL
˙ − (T¯(1) − T (1) )vdx− 0
h
(1)
+ b2 T1 v˙
iL 0
RL
(2)
RL
− b2 T1 vdx ˙ − (T¯ − T )ϕ˙ = 0. 0
0
(37) Applying the standard localization procedure, the following set of ordinary differential equations (ODEs) is obtained: u˙ 0 : v˙ : ϕ˙ :
(1) (1) (2) P¯ (1) − PM − PT + c2 P1 = 0, (2) −T¯(1) + T (1) − b2 T1 = 0, (1) (1) (2) (1) ¯ (1) − MM M + T¯ − T − MT + c2 M1 + b2 T1 = 0,
(38)
with kinematic boundary conditions for {u0 , v, ϕ} given at {0, L} or alterna12
ACCEPTED MANUSCRIPT
tively via static conditions: (1) prescribe u0 or P¯ − PM − PT + c2 P1 = 0, (1) prescribe u0 or − c2 P1 = 0, (1) prescribe v or T¯ − T + b2 T1 = 0, prescribe v (1) or − b2 T1 = 0, ¯ − MM − MT + c2 M1(1) + b2 T1 = 0, prescribe ϕ or M prescribe ϕ(1) or − c2 M1 = 0.
(39)
¯ (1) − T¯, distributed transversal Note that distributed couples are m = −M load py = −T¯(1) and distributed longitudinal load px = −P¯ (1) , see [33] for a detailed elaboration. The transversal and the longitudinal forces in a crosssection correspond to T¯ and P¯ , respectively, while the couple in a cross¯. section corresponds to M Solution of the above ODE system can be now performed in two ways. First one is a standard one and is most frequently used, see [33] for an example. The second one, as proposed in this paper leads to a more simpler system and will be favoured in the example section. Standard approach: Consider just introduced external loads, Eqs. (31,32) and Eqs. (34). The governing ODEs Eqs. (38) can be transformed by straightforward manipulation into the final form: u˙ 0 : v˙ : ϕ˙ :
(4)
(2)
(1)
−c2 ϕ(4) kES + ϕ(2) kES + c2 u0 kEA − u0 kEA − PT − px = 0, kGA b2 (v (4) − ϕ(3) ) − kGA (v (2) − ϕ(1) ) − py = 0, (4) c2 (ϕ(4) kEI − u0 kES ) + b2 kGA (−ϕ(1) + v (2) ) + kGA (ϕ − v (1) )− (2) (1) −ϕ(2) kEI + u0 kES − MT − m = 0, (40)
13
ACCEPTED MANUSCRIPT
while the accompanying boundary conditions are: (1) (3) prescribe u0 or − c2 (ϕ(3) kES − u0 kEA ) + ϕ(1) kES − u0 kEA − PT + P¯ = 0, (1) (2) prescribe u0 or ϕ(2) kES − u0 kEA = 0, prescribe v or − b2 kGA (ϕ(2) − v (3) ) + kGA ϕ − kGA v (1) + T¯ = 0, prescribe v (1) or b2 kGA (ϕ(1) − v (2) ) = 0, (3) prescribe ϕ or c2 (ϕ(3) kEI − u0 kES ) − b2 kGA (ϕ(1) − v (2) )− (1) ¯ = 0, −ϕ(1) kEI + u0 kES − MT + M (2) (1) 2 (2) prescribe ϕ or c (−ϕ kEI + u0 kES ) = 0. (41) It can be easily seen, the resulting system of ODE is significantly coupled. Alternative approach: Alternatively to the usual practice, we do not transform the problem to the form involving displacements and the angle as unknowns. Instead, we rather choose to work with forces and moments as unknowns. First, note that from Eqs. (34) follows: (1)
M1 = MM ,
T1 = T (1) ,
(1)
P1 = PM .
(42)
Applying the above result, the Eqs. (38) take the following form u˙ 0 : v˙ : ϕ˙ :
(3) (1) (1) P¯ (1) − PM − PT + c2 PM = 0, −T¯(1) + T (1) − b2 T (3) = 0, (3) (1) (1) ¯ (1) − MM + T¯ − T − MT + c2 MM + b2 T (2) = 0, M
(43)
with boundary conditions : (2) prescribe u0 or P¯ − PM − PT + c2 PM = 0, (1) (1) prescribe u0 or − c2 PM = 0, prescribe v or T¯ − T + b2 T (2) = 0, prescribe v (1) or − b2 T (1) = 0, ¯ − MM − MT + c2 M (2) + b2 T (1) = 0, prescribe ϕ or M M (1) (1) 2 prescribe ϕ or − c MM = 0.
(44)
Obviously, Eqs. (43)1,2 are decoupled and can be easily solved. Subsequently, with the shear and axial forces defined, the solution of Eq. (43)3 is straightforward.
14
ACCEPTED MANUSCRIPT
Some of unknown integration constants follow from the static boundary conditions. This is not always the case. For example, if the doubly clamped beam is analysed, not all integration constants can be determined from static boundary conditions. Unknown derivatives of displacement functions can be now easily obtained from the solution of the linear system Eqs. (34). Application of the kinematic boundary conditions is used to determine integration constants in this operation as well as the remaining integration constants (if any) from the solving of static quantities. From another point of view, the problem can be also imagined as the set of 6 ODE, Eqs. (43) and Eqs. (34)1,3,5 , six unknown functions {PM , T, MM , u0 , v, ϕ}, accompanied by the suitable set of boundary conditions Eqs. (44). 4. Examples 4.1. Functionally graded cantilever beam loaded by the non-homogeneous temperature field Let us consider a cantilever nanobeam of length L and of the rectangular cross-section HxB. The longitudinal axis is denoted as x while bending is assumed to take place in x − y plane. Young’s and shear modulus are linear functions of the coordinate y: EFG = E0 + E1 y,
GFG = G0 + G1 y,
(Eα)FG = Eα0 + Eα1 y,
(45)
where E0 , E1 , G0 , G1 are some material constants. The product (Eα)FG also assumes a linear variation along the coordinate y defined by additional material parameters Eα0 , Eα1 . The temperature field is non-homogeneous and is defined by the quadratic distribution θ = θ(y) along the height of the nanobeam: ∆θ(y) = t0 + t1 y + t2 y 2 . (46) Additionally, a case θ = θ(x, y) where the above distribution is scaled by the longitudinal coordinate x is considered as well: ∆θ(x, y) = (t0 + t1 y + t2 y 2 )x.
(47)
The beam is not mechanically loaded in any way and is clamped at the right ¯ (0) = 0 and the Eqs. (43) take the end. Consequently,P¯ (0) = T¯(0) = 0, M
15
ACCEPTED MANUSCRIPT
following form (1)
(1)
(3)
−PM − PT + c2 PM = 0, T (1) − b2 T (3) = 0, (1) (1) (3) −MM + T¯ − T − MT + c2 MM + b2 T (2) = 0.
(48)
Boundary conditions, Eqs. (39), are: at x = 0 (2) −PM − PT + c2 PM = 0, −T + b2 T (2) = 0, (2) −MM − MT + c2 MM + b2 T (2) = 0, (1) (1) −c2 PM = 0, −b2 T (1) = 0, −c2 MM = 0 at x = L u0 = 0, v = 0, ϕ = 0, (1) (1) −c2 PM = 0, −b2 T (1) = 0, −c2 MM = 0.
(49)
Case θ = θ(y) : In the case when temperature is defined by the Eq. (46), solutions of Eq. (48) are: PM = cex/c c1 − ce−x/c c2 + c3 , T = c4 + 21 be−x/b (−1 + e2x/b )c5 + 12 b2 e−x/b (−1 + ex/b )2 c6 , x x MM = c7 + 12 ce− c (−1 + e2x/c )c8 + 12 c2 e− c (−1 + ex/c )2 c9 + x 2x x 2x + 12 e− c (−c + ce c − 2xex/c )c4 − 12 b2 e− c (−c + ce c − 2xex/c )c6 .
(50)
Solving for the integration constants, it is obtained: c1 = c2 = c4 = c5 = c6 = c8 = c9 = 0, 1 c3 = BHEα0 t0 + 12 BH 3 (Eα1 t1 + Eα0 t2 ), 1 1 3 c7 = − 12 BH (Eα0 t1 + Eα1 t0 ) − 80 BH 5 Eα1 t2 ,
(51)
thus collapsing the solution to: 1 PM = c3 = BHEα0 t0 + 12 BH 3 (Eα1 t1 + Eα0 t2 ), T = 0, 1 1 MM = c7 = − 12 BH 3 (Eα0 t1 + Eα1 t0 ) − 80 BH 5 Eα1 t2 ,
(52)
i.e, both PM and MM are independent of the coordinates x while the shear 16
ACCEPTED MANUSCRIPT
force vanishes. Now, from Eqs. (34) it is: (1)
PM = −ϕ(1) kES + u0 kEA , 0 = (v (1) − ϕ)kGA , (1) MM = ϕ(1) kEI − u0 kES ,
(53)
along with boundary conditions u0 (L) = 0, v(L) = 0, ϕ(L) = 0. Integrating for the unknown functions, it is obtained: M E0 +MM E1 )x + c10 , u0 = − 12(P BH(E 2 H 2 −12E 2 ) 1
0
2
2
M E0 +PM E1 H )x v = − (6(12M + c11 + xc12 , BH 3 (−12E 2 +E 2 H 2 )) 0
ϕ=
1
2 M E0 +PM E1 H )x − (12(12M BH 3 (−12E02 +E12 H 2 ))
(54)
+ c12 ,
with integration constants: c10 =
2 2 2 (12(PM E0 +MM E1 )L) M E0 L +6PM E1 H L , c11 = 72M12BE , 2 H 3 −BE 2 h5 BH(−12E02 +E12 H 2 ) 0 1 2 M E0 +PM E1 H )L c12 = 12(12M . BH 3 (−12E02 +E12 H 2 )
(55)
Thus, the axial displacement and the bending angle are linearly distributed, while the transverse displacement takes the parabolic form. Case θ = θ(x, y) : In the case when temperature is defined by the Eq. (47), solutions of Eqs. (48) are: 1 B(12Eα0 Ht0 + Eα1 H 3 t1 + Eα0 H 3 t2 )x + cex/c c1 − ce−x/c c2 + c3 , PM = 12 T = c7 + 21 be−x/b (−1 + e2x/b )c8 + 21 b2 e−x/b (−1 + ex/b )2 c9 , 1 MM = 960 cBH 3 (−1 + ex/c )2 (1 − e−2x/c )(20Eα1 t0 + 20Eα0 t1 + 3Eα1 H 2 t2 )− 1 cBH 3 (1 + e−2x/c )(−1 + e2x/c )(20Eα1 t0 + 20Eα0 t1 + 3Eα1 H 2 t2 )+ 960 1 BH 3 (20Eα1 t0 + 20Eα0 t1 + 3Eα1 H 2 t2 )(−ce−x/c + cex/c − 2x) + c4 + 480 1 −x/c ce (−1 + e2x/c )c5 + 21 c2 e−x/c (−1 + ex/c )2 c6 + 21 e−x/c (−c + ce2x/c − 2ex/c x)c7 − 2 1 2 −x/c be (−c + ce2x/c − 2ex/c x)c9 , 2 (56)
17
ACCEPTED MANUSCRIPT
Solving for the integration constants, it is obtained: 2
α1 H t1 +Eα0 H c1 = − BH(12Eα0 t0+E 12(1+eL/c )
c2 c3 c4 c5 c6 c7
2t ) 2
,
BH(12Eα0 t0 +Eα1 H 2 t1 +Eα0 H 2 t2 ) 1 B(12Eα0 Ht0 + Eα1 H 3 t1 − 12 12(1+eL/c ) 1 (12BEα0 ht0 x + BEα1 H 3 t1 x + BEα0 H 3 t2 x), 12 3 2 L/c )+x+eL/c x)) α1 H t2 )(−c+ce − BH (20Eα1 t0 +20Eα0 t1 +3E , 240(1+eL/c ) 1 3 2 BH (20Eα1 t0 + 20Eα0 t1 + 3Eα1 H t2 ), 240 BH 3 (−1+eL/c )(20Eα1 t0 +20Eα0 t1 +3Eα1 H 2 t2 ) − , 240(1+eL/c )
= = = = = = c8 = c9 = 0,
+ Eα0 H 3 t2 ),
(57) thus collapsing the solution Eqs. (48) to: −x/c
2
2
L/c
2x/c
x/c
L/c+x/c
H t2 )(−ce +ce −2e x−2xe ) PM = − BHe (12Eα0 t0 +Eα1 H t1 +Eα012(1+e L/c ) T = 0, 3 −x/c 2 t )(−ceL/c +ce2x/c −2xex/c −2xeL/c+x/c ) 2 MM = BH e (20Ea1t0 +20Ea0t1 +3Ea1H . 240(1+eL/c )
(58) Differently from the previous case, PM and MM are now dependent on the coordinate x, while the shear force vanishes once more. Again, from Eqs. (34) and with boundary conditions u0 (L) = 0, v(L) = 0, ϕ(L) = 0,it is obtained: L−x c +c2 (−ex/c )+x2 (eL/c +1))(20E0 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 H 2 t1 )−E1 H 2 (20Eα0 t1 +3Eα1 H 2 t2 +20Ea1t0 )) 20(eL/c +1)(12E02 −E12 H 2 ) L−x 2 3 L/c 2 2 x/c c (−c (c+x)e −c e (x−c)+ 3 x (e +1))5E1 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 H 2 t1 )−3E0 (20Eα0 t1 +3Eα1 H 2 t2 +20Eα1 t0 ) 5(eL/c +1)(12E02 −E12 H 2 ) L−x 2 2 x/c 2 L/c 2 x(−c e c −c e +x (e +1))(5E1 (Eα0 H t2 +12Eα0 t0 +Eα1 H 2 t1 )−3E0 (20Eα0 t1 +3Eα1 H 2 t2 +20Eα t0 )) 12 11 5(eL/c +1)(12E02 −E12 H 2 ) L−x 2 2 x/c 2 L/c 2 2 2 c (−c e +c (−e )+x (e +1))(5E1 (Eα0 H t2 +12Eα0 t0 +Eα1 H t1 )−3E0 (20Eα0 t1 +3Eα1 H t2 +20Eα1 t0 )) 12 5(eL/c +1)(12E02 −E12 H 2 )
u0 =
(−c2 e
+ c10 ,
v=
+
+c x+c ,
ϕ=
+c ,
(59) with integration constants: c10 =
(c−L)(c+L)(20E0 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 H 2 t1 )−E1 H 2 (20Eα0 t1 +3Eα1 1H 2 t2 +20Eα1 t0 )) , 240E02 −20E12 H 2 −L/2c
L/2c
c11 =
(3c3 eL/2c −e−L/2c −3c2 L+2L3 )(5E1 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 H 2 t1 )−3E0 (20Eα0 t1 +3Eα1 H 2 t2 +20Eα1 t0 )) e
c12 =
+e
180E02 −15E12 H 2 (c−L)(c+L)(5E1 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 a1H 2 t1 )−3E0 (20Eα0 t1 +3Eα1 H 2 t2 +20Eα1 t0 )) 60E02 −5E12 H 2
(60) Therefore, if these two cases with different temperature distributions are
18
.
,
ACCEPTED MANUSCRIPT
compared, it can be said that in the beam that is not loaded by mechanical forces, shear forces do not occur. As expected, longitudinal temperature variation leads toward more complex solutions. In particular, the axial displacement in the case of θ = θ(y) is distributed linearly, while in the second case significantly more complex distribution involving exponential functions take place. 4.2. Simply supported functionally graded beam loaded by non-homogeneous temperature field and mechanical quantities Let us consider a simply supported nanobeam of L = 100 nm and of the rectangular cross-section H = 20 nm, B = 10 nm. As in the previous example, the beam longitudinal axis is denoted with x and deformation is assumed to take place in x − y plane. At the left end, the beam elongation is constrained in all directions, ie. u(0) = v(0) = 0, while at the right end beam dilatation along x axis is free, so v(L) = 0. Shear correction factor is χ = 5/6. The nanobeam is made of Poly(methyl methacrylate) (PMMA) matrix. Young’s modulus of PMMA is EPMMA = 2.5 GPa, shear modulus GPMMA = 0.93 GPa, specific density ρPMMA = 1190 kg/m3 and the coefficient of thermal expansion αPMMA = 70 · 10−6 K−1 . As a reinforcement, single walled carbon nanotubes (CNT) are used. CNTs Young’s modulus is ECNT = 1 TPa, GCNT = 0.5 TPa and αCNT = −1.2 · 10−5 K−1 . Note that the negative sign of CTE describes contraction of the CNT with the temperature increase. Nanotubes are dispersed with approximately linear variation along CNT the beam height, as depicted in Fig. 1, i.e. 30 CNTs in total. Therefore, CNT volume fraction is: PMMA
CNT PMMA
Figure 1: The cross section of the simply supported nanobeam
VCNT = 2
|y| max , V H CNT
19
(61)
ACCEPTED MANUSCRIPT
max where VCNT represent maximal volume fraction of CNTs at the y = ±H/2. max In this case VCNT = 0.1. Consequently, according to the Eq. (1), PMMA volume fraction must be:
VPMMA = 2
H − |y| Vmax . H
(62)
Three loading cases are considered: Case a. The beam is thermally loaded by the homogeneous temperature field, ∆θ = aa = 30 K. Case b. The beam is thermally loaded by the linear variation of temperature field along the longitudinal axis, ∆θ = ab x = 0.3x K. Case c. The beam is thermally loaded by the linear variation of temperature field along the longitudinal axis, ∆θ = ac = 0.3x K and by the downward distributed load q = 5 nN/nm. According to the above defined loading, it is px = 0, m = 0. In Cases a and b, py = 0, while in the Case c it is py = −q. Eventually, Eqs. (43) collapse to: u˙ 0 : v˙ : ϕ˙ :
(1) (3) P¯ (1) − PM + c2 PM = 0, −py + T (1) − b2 T (3) = 0, (1) (3) −MM + T¯ − T + c2 MM + b2 T (2) = 0.
(63)
¯ (0) = M ¯ (L) = 0 and due Since the beam is simply supported, P¯ (L) = 0, M to boundary conditions, Eqs. (39), it is : at x = 0 u0 = 0, v = 0 (2) −MM − MT + c2 MM + b2 T (2) = 0, (1) (1) −c2 PM = 0, −b2 T (1) = 0, −c2 MM = 0, at x = L (2) −PM − PT + c2 PM = 0, v = 0, (1) −MM − MT + c2 MM + b2 T (1) = 0, (1) (1) −c2 PM = 0, −b2 T (1) = 0, −c2 MM = 0. 20
(64)
ACCEPTED MANUSCRIPT
Upon determination of PM , T, MM , the displacements and the rotation can be evaluated. Temperature distribution does not depend on the longitudinal coordinate, so ∂x kθA = 0, ∂x kθS = 0. Consequently, system of equations provided by Eqs. (34) collapses to: (1)
MM = ϕ(1) kEI − u0 kES , T = v (1) kGA − ϕkGA = (v (1) − ϕ)kGA , (1) PM = −ϕ(1) kES + u0 kEA .
(65)
With the above equations at hand, for each of the three described loading cases, the solutions are obtained as follows. 4.2.1. Case a. The thermal axial force and the thermal bending moment, according to the Eqs. (31)4 ,(32)7 , are PT =
R
σT da = − (Eα)FG (θ − θ0 )da = 45.6 nN, R
A
A
MT = − σT yda = (Eα)FG (θ − θ0 )yda = 8.1 nNnm. R
R
A
A
(66)
With the above values, solution of the ODE system Eqs. (70) is: PM = −PT = −45.6 nN, T = 0, MM = −MT = −8.1 nNnm.
(67)
Consequently, both PM and MM do not change along the beam axis, while shear forces do not occur. Now, displacements and rotation fields can be evaluated from Eqs. (65): ES MT x = −23.98 · 10−5 x, u = − KKEIEAPTK+K 2 EI −K ES
v= ϕ=
(KEA LMT +KES LPT )x−(KEA MT +KES PT )x2 = 1.14 · 10−6 x 2 2(KEA KEI −KES ) (L−2x)(KEA MT +KEs PT ) = 1.14 · 10−6 − 2.272 · 10−8 x. 2 2(KEA KEI −KES )
− 1.136 · 10−8 x2 ,
(68) These results are presented in Figs. 2-4. Evidently, nonlocal effects play no role in the case of homogeneous temperature field. Please note that BernoulliEuler and Timoshenko formulations give identical results for a such case. Negative sign in the axial displacement field indicate that the effect of thermal 21
ACCEPTED MANUSCRIPT
contraction of carbon nanotubes dominates - leading toward contraction of the nanobeam. 0.000 -0.005 uHnmL
-0.010 -0.015 -0.020 -0.025
0
20
40
60
L HnmL
80
100
Figure 2: Longitudinal displacements of the simply supported nanobeam in the case of the constant temperature field
0.00003 0.000025 0.00002 vHnmL 0.000015 0.00001 5. ´ 10 -6 0
0
20
40
60
L HnmL
80
100
Figure 3: Transversal displacements of the simply supported nanobeam in the case of the constant temperature field
4.2.2. Case b. In this case, the thermal axial force and the thermal bending moment Eqs. (31)4 ,(32)7 , are not constant along the nanobeam: PT =
R
σT da = − (Eα)FG (θ − θ0 )da = PT0 x, R
A
A
MT = − σT yda = (Eα)FG (θ − θ0 )yda = MT0 x, R
R
A
A
22
(69)
ACCEPTED MANUSCRIPT
1.5 ´ 10 -6 1. ´ 10 -6 5. ´ 10 -7
j HradL
0 -5. ´ 10 -7 -1. ´ 10 -6 -1.5 ´ 10 -6
0
20
40
60
80
L HnmL
100
Figure 4: Rotation angle of the simply supported nanobeam in the case of the constant temperature field
where
max max b PT0 = − BHa + 2E2 α2 ), − E2 α2 VCNT (E1 α1 VCNT 2 2 max BH VCNT ab MT0 = (E1 α1 − E2 α2 ). 12
(70)
Consequently, both PT and MT along the beam axis vary in the linear manner. With the above values, solution of the ODE system Eqs. (70) is:
L
L
x
x
PT0 −xe c −ce c − c +ce c −x
PM = MM =
MT0
L e c +1 L L x x −xe c −ce c − c +ce c −x
(71) ,
L
e c +1
T = 0. Once more, shear forces do not occur. Having in mind that KEA , KES and KEI are not functions of the longitudinal coordinate x, displacements and rotation fields can be evaluated from Eqs. (71):
u=
L
L
x
x
L
2 e c +1
h
v=
L
2c2 e c −2c2 e c − c +2c2 −2c2 e c +x2 e c +x2 (kEI PT0 +kES MT0 )
L
x
2 (kEA kEI −kES )
L
L
L
6L e c +1
ϕ=
i
L
6c3 (Le c −L+Le c −2xe c −LeL/c−x/c +2x)+(L3 xe c +L3 x−Lx3 e c −Lx3 ) (kEA MT0 +kES PT0 )
L
x
L
x
2 (kEA kEI −kES )
L
L
−12c3 e c +12c3 +6c2 Le c +6c2 Le c − c +L3 e c +L3 −3Lx2 e c −3Lx2 (kEA PM0 +kES PT0 )
L
6L e c +1
2 (kEA kEI −kES )
. (72)
23
ACCEPTED MANUSCRIPT
0 -10 Pm HnNL
-20 -30
c=1
-40
c=20
-50
c=1000
-60 0
20
40
60
80
100
LHnmL
Figure 5: Variation of Pm along the simply supported nanobeam under the linear temperature field for increasing values of nonlocal parameter c 0
Mm HnNnmL
-2 -4 c=1
-6
c=20
-8
c=1000
-10 -12
0
20
40
60
80
100
LHnmL
Figure 6: Variation of Mm along the simply supported nanobeam under the linear temperature field for increasing values of nonlocal parameter c 0.000
u HnmL
-0.005 c=1
-0.010
c=20 c=1000
-0.015 0
20
40
60
80
100
LHnmL
Figure 7: Longitudinal displacements of the simply supported nanobeam under the linear temperature field for increasing values of the nonlocal parameter c
24
ACCEPTED MANUSCRIPT
0.00002
vHnmL
0.000015 c=1
0.00001
c=20
5. ´ 10-6
c=1000
0
0
20
40
60
80
100
LHnmL
Figure 8: Transversal displacements of the simply supported nanobeam under the linear temperature field for increasing values of the nonlocal parameter c
j HradL
5. ´ 10-7 0
c=1 c=20
-5. ´ 10-7
c=1000
-1. ´ 10-6 0
20
40
60
80
100
LHnmL
Figure 9: Rotation angle of the simply supported nanobeam under the linear temperature field for increasing values of the nonlocal parameter c
25
ACCEPTED MANUSCRIPT
Note that PM and MM are described by the same class of curves and differ by a factor PT0 /MT0 , compare Figs. 5,6. Nonlocal influence on the transverse displacement v is manifested through a slight increase of the maximal displacement and shift of the extremum away from from the midpoint span. 4.2.3. Case c. Following the same procedure as before, the solution can be readily obtained. The thermal axial force and the thermal bending moment take same values as in the Case b, PT = PT0 x, MT = MT0 x. Now, the solution of the ODE system, Eqs. (70) is:
PM =
x
L
L
x
PT0 −xe c −ce c − c +ce c −x
T =− MM =
py
,
L e c +1 L x L L L−x L2 e b +L2 +2bLe b −2Lx(e b +1)−4be b −2bLe b b +4b
e 2L(
L/b +1
2L
1
L e b +1
L − xc c
2L e c −1
L
L
,
) L
x
L
L
L
2L
L
py e c L2 c(e b + e c + e b + c + 1)
L
L
L
2L
+e L2 c(e b + e b + c + e c + 1) + 2Lc2 (e b − e c − e b + c + 1) L L L L x L L L x L L +4bc(e c − c (e b − e b + c + e c − 1) − e c (e b + e c − e b + c − 1)) L
+2b(2x − L)(e b + e
L
2L c
L
2L c
−eb+
L
L
x
L
− 1) + xL(x − L)(e b − e L
x
L
x
2L c
L
−eb+
2L c
+ 1)
x
+2MT0 L c(1 − e c )(e c − c + e b + c − c − e b + c − e c ) L
+x(e b − e
2L c
L
−eb+
2L c
+ 1)
. (73)
Distribution of PM is identical to the result obtained in Case b. Unfortunately, resulting functions that describe displacements and the rotation take rather lengthy form and due to economy in space, the formal results are not presented. Instead, these are presented in the graphical form, thus illustrating main points. As visible from Fig. 10 and Eqs. (73)2 , distribution of T is insensitive of the nonlocal parameter c. From Fig. 11 it can be noticed that the influence of the nonlocal parameter b is much more pronounced than the nonlocal parameter c. The bending moment MM is pretty much insensitive of the nonocal parameter b. The axial displacement at the midpoint of the nanobeam, Fig. 12, does not depend on b too, while transversal displacements and rota-
26
ACCEPTED MANUSCRIPT
tion angle, Figs. 13 and 14, depend on both nonlocal parameters. Sensitivity with respect to both b and c is approximately equal when the transversal displacement is considered, while influence on nonlocal parameters on the rotation angle manifest different trends. 150
b
100
50
0 4 2 T 0
50 100
0
c
150
Figure 10: T at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
0
50 c 100 150 6000 5500 5000 M M 4500 0 50
100
b
150
Figure 11: MM at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
5. Conclusions In the paper at hand a model for nanoscale thermomechanical analysis of beams was presented. In order to properly describe effects occurring at small scales, the extended version of Erignen’s original nonlocal formulation was proposed. The nanobeam model was based on Timoshenko formulation. 27
ACCEPTED MANUSCRIPT
- 0.011 - 0.012 - 0.013 - 0.014 u
0 50 c
0 50
100
100 b 150 150 0
50 c
100 150
150
100 b
50
0
- 0.011 - 0.012 - 0.013 - 0.014 u
Figure 12: The axial displacement at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
28
ACCEPTED MANUSCRIPT
- 0.9 v 1.0 - 1.1 150
150 100
100 c
50
50
b
00
- 0.9 - 1.0
150
- 1.1 150
100 c
v
100
50
50
b
0 0
Figure 13: The transversal displacement at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
29
ACCEPTED MANUSCRIPT
c
0
50
100 150 0.0002 j 0.0001 0.0000 150
100
50
b c 150
100
50
0
0
0.0002 0.0001 0.0000 j 150
100 b
50
0
Figure 14: The rotation angle at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
30
ACCEPTED MANUSCRIPT
The approach was verified on two examples, thus establishing two new benchmarks in the field. The main findings point out necessity of employment of the Timoshenko formulation instead of the Bernoulli-Euler formulation if the loading, either mechanical or thermal, takes other form than the most simplest. Thus, although the Bernoulli-Euler formulation or isothermal models are simpler to implement, at the same time they offer reduced possibilities of mechanical structural modelling at the nano level. In that way, the present formulation surpasses previous findings providing more general solutions and wider field of application. Moreover, the strict thermodynamical setting assures that second law of thermodynamics is fulfilled for all possible scenarios, hence providing the physically sound formulation. 6. Acknowledgements This work has been fully supported by Croatian Science Foundation under the project no. 6876 - Assessment of structural behaviour in limit state operating conditions. This support is gratefully acknowledged. References References [1] Liew K, Lei Z, Zhang L. Mechanical analysis of functionally graded carbon nanotube reinforced composites: A review. Composite Structures 2015;120(0):90 –7. [2] Kiani K. Longitudinal and transverse instabilities of moving nanoscale beam-like structures made of functionally graded materials. Compos Struct 2014;107:610–9. [3] Nazemnezhad R, Hosseini-Hashemi S. Nonlocal nonlinear free vibration of functionally graded nanobeams. Compos Struct 2014;110(1):192–9. [4] Hosseini-Hashemi S, Nazemnezhad R, Bedroud M. Surface effects on nonlinear free vibration of functionally graded nanobeams using nonlocal elasticity. Appl Math Modell 2014;38(14):3538–53. [5] Ke LL, Yang J, Kitipornchai S. Nonlinear free vibration of functionally graded carbon nanotube-reinforced composite beams. Compos Struct 2010;92(3):676–83. 31
ACCEPTED MANUSCRIPT
[6] Simsek M. Non-linear vibration analysis of a functionally graded Timoshenko beam under action of a moving harmonic load. Compos Struct 2010;92(10):2532–46. [7] Yas M, Heshmati M. Dynamic analysis of functionally graded nanocomposite beams reinforced by randomly oriented carbon nanotube under the action of moving load. Appl Math Modell 2012;36(4):1371–94. [8] Simsek M, Yurtcu H. Analytical solutions for bending and buckling of functionally graded nanobeams based on the nonlocal Timoshenko beam theory. Compos Struct 2013;97:378–86. [9] Eltaher M, Khairy A, Sadoun A, Omar FA. Static and buckling analysis of functionally graded Timoshenko nanobeams. Appl Math Comput 2014;229:283–95. [10] Eltaher M, Emam S, Mahmoud F. Static and stability analysis of nonlocal functionally graded nanobeams. Compos Struct 2013;96:82–8. [11] Barretta R, Feo L, Luciano R, Marotti de Sciarra F. Variational formulations for functionally graded nonlocal Bernoulli-Euler nanobeams. Compos Struct 2015;129:80–9. [12] Marotti de Sciarra F, Canadija M, Barretta R. A gradient model for torsion of nanobeams. C R Mec 2015;343(4):289–300. [13] Barretta R, Feo L, Luciano R. Torsion of functionally graded nonlocal viscoelastic circular nanobeams. Composites Part B: Engineering 2015;72:217–22. [14] Erisken C, Kalyon D, Wang H. Functionally graded electrospun polycaprolactone and β -tricalcium phosphate nanocomposites for tissue engineering applications. Biomaterials 2008;29(30):4065–73. [15] Filiz S, Aydogdu M. Wave propagation analysis of embedded (coupled) functionally graded nanotubes conveying fluid. Compos Struct 2015;132:1260–73. [16] Friedrich K, Zhang Z, Schlarb A. Effects of various fillers on the sliding wear of polymer composites. Compos Sci Technol 2005;65(15-16 SPEC. ISS.):2329–43. 32
ACCEPTED MANUSCRIPT
[17] Kwon H, Bradbury C, Leparoux M. Fabrication of functionally graded carbon nanotube-reinforced aluminum matrix composite. Adv Eng Mater 2011;13(4):325–9. [18] Canadija M, Barretta R, de Sciarra M. A gradient elasticity model of Bernoulli-Euler nanobeams in non-isothermal environments. submitted ????;. [19] Shen HS. Postbuckling of nanotube-reinforced composite cylindrical shells in thermal environments, part I: Axially-loaded shells. Compos Struct 2011;93(8):2096–108. [20] Shen HS, Xiang Y. Nonlinear vibration of nanotube-reinforced composite cylindrical shells in thermal environments. Comput Meth Appl Mech Eng 2012;213-216:196–205. [21] Khater M, Eltaher M, Abdel-Rahman E, Yavuz M. Surface and thermal load effects on the buckling of curved nanowires. Engineering Science and Technology, an International Journal 2014;17(4):279 –83. [22] Zhang YW, Chen J, Zeng W, Teng YY, Fang B, Zang J. Surface and thermal effects of the flexural wave propagation of piezoelectric functionally graded nanobeam using nonlocal elasticity. Comput Mater Sci 2015;97:222–6. [23] Zenkour A, Abouelregal A. Nonlocal thermoelastic nanobeam subjected to a sinusoidal pulse heating and temperature-dependent physical properties. Microsyst Technol 2014;:1–10. [24] Zenkour A, Abouelregal A, Alnefaie K, Abu-Hamdeh N, Aifantis E. A refined nonlocal thermoelasticity theory for the vibration of nanobeams induced by ramp-type heating. Appl Math Comput 2014;248:169–83. [25] Alshorbagy A. Temperature effects on the vibration characteristics of a functionally graded thick beam. Ain Shams Engineering Journal 2013;4(3):455–64. [26] JafarSadeghi-Pournaki I, Zamanzadeh M, Madinei H, Rezazadeh G. Static pull-in analysis of capacitive FGM nanocantilevers subjected to thermal moment using Eringen’s nonlocal elasticity. Int J Eng Trans A 2014;27(4):633–42. 33
ACCEPTED MANUSCRIPT
[27] Bafekrpour E, Simon G, Habsuda J, Naebe M, Yang C, Fox B. Fabrication and characterization of functionally graded synthetic graphite/phenolic nanocomposites. Mater Sci Eng, A 2012;545:123–31. [28] Bafekrpour E, Simon G, Yang C, Chipara M, Habsuda J, Fox B. Functionally graded carbon nanofiber-phenolic nanocomposites for sudden temperature change applications. Polymer (United Kingdom) 2013;54(15):3940–8. [29] Mehrabadi S, Aragh B, Khoshkhahesh V, Taherpour A. Mechanical buckling of nanocomposite rectangular plate reinforced by aligned and straight single-walled carbon nanotubes. Composites Part B: Engineering 2012;43(4):2031–40. [30] Shen HS, Noda N. Postbuckling of pressure-loaded fgm hybrid cylindrical shells in thermal environments. Compos Struct 2007;77(4):546–60. [31] Ahankari S, Kar K. Processing of styrene butadiene rubber-carbon black nanocomposites with gradation of crosslink density: Static and dynamic mechanical characterization. Mater Sci Eng, A 2008;491(1-2):454–60. [32] Eringen A. On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves. J Appl Phys 1983;54(9):4703–10. [33] Marotti De Sciarra F, Barretta R. A gradient model for Timoshenko nanobeams. Physica E 2014;62:1–9. [34] Eringen A. Nonlocal Continuum Field Theories. Springer Verlag; 2002. [35] Coleman B, Gurtin M. Thermodynamics with internal state variables. J Chem Phys 1967;47:597–613.
34
S0263-8223(15)00873-9 http://dx.doi.org/10.1016/j.compstruct.2015.09.030 COST 6875
To appear in:
Composite Structures
Please cite this article as: Čanađija, M., Barretta, R., de Sciarra, F.M., On Functionally Graded Timoshenko Nonisothermal Nanobeams, Composite Structures (2015), doi: http://dx.doi.org/10.1016/j.compstruct.2015.09.030
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
On Functionally Graded Timoshenko Nonisothermal Nanobeams Marko Čanađijaa,∗, Raffaele Barrettab , Francesco Marotti de Sciarrab a
Faculty of Engineering,Department of Engineering Mechanics,University of Rijeka,Vukovarska 58,51000 Rijeka,Croatia b Department of Structures for Engineering and Architecture,University of Naples Federico II,Via Claudio 21,80121 Naples,Italy
Abstract The paper deals with the mechanics of nanobeams within nonisothermal environments. Mechanical behaviour is assumed to be time-independent and suitable for the statical beam deformation. The nanobeam mechanics is based on the Timoshenko kinematics, augmented by nonlocality effects as advocated by Eringen. Following a thermodynamic approach, a sound framework oriented towards anisotropic materials is developed. Such nanobeam model is especially suited for functionally graded materials. The proposed procedure is tested on two examples. The first example investigates mechanical behaviour in the more general way for the different temperature distributions, while the second one is oriented toward specific carbon nanotube Poly(methyl methacrylate) nanocomposite. Keywords: Timoshenko nanobeams, nonlocal thermoelasticity, carbon nanotubes 1. Introduction A great interest of researchers is nowadays directed toward carbon nanotubes (CNT) reinforced polymers. This is of course related to the superior CNT mechanical properties, but additional benefits can be further achieved by the functionally grading of the composite [1]. The theoretical research of the subject can be divided into three main directions - static and buckling, ∗
Corresponding author. Tel.: +385-51-651-496; Fax.: +385-51-651-490 Email address: [email protected] (Marko Čanađija )
Preprint submitted to Composite Structures
September 23, 2015
ACCEPTED MANUSCRIPT
dynamic and vibrations and nonlinear analysis. Remarkably, only a tiny fraction of published papers addresses nonisothermal behaviour. Within the isothermal context, previous non-local functionally graded (FG) beam models were developed to address vibration issues [2, 3, 4, 5, 6, 7]. Interests in the vibrational response of a nanobeam lies in the frequent application in nanosenors. In particular, the resonant frequency is exploited in mass sensing at the nanoscale. Another class of sensors are based on the classical beam statics. Hence, bending and buckling of a non-local FG nanobeam was investigated in [8, 9] for Timoshenko formulation and Bernoulli-Euler nanobeam in [10, 11]. Torsion was investigated in [12]. Extension to viscoelastic behaviour was recently proposed in [13]. Beside nanosensorics, application of FG nanocomposites is very wide; it seems that tissue engineering can significantly benefit from such materials, resembling composition and structure of the bone tissue much better than the conventional materials, [14]. FG nanotubes can be used for fluid transport in the nano electromechanical systems [15]. Different grading distributions can be conveniently obtained by the centrifugation techniques [16], resulting in polymer composites with superior wear resistance. Another frequent manufacturing process of FG nanocomposites uses assembling of different layers of materials. For example, aluminium - CNT composite is produced from layers containing different fractions of CNT [17]. In this particular case, Vickers hardness of each particular layer was measured with results varying by a factor of 7. However, when influence of temperature on the nanomechanics is considered, there are only a few papers available. Bernoulli-Euler kinematics of nonisothermal nanobeams for isotropic materials is analyzed in [18]. Previous FG material thermomechanical applications at nanoscale involve bifurcation analysis of functionally graded (FG) shells [19] and nonlinear vibrations of FG shells [20]. Thermal buckling of curved nanowires, yet with a homogeneous temperature field was investigated in [21]. Nonlocality, by means of surface effects were accounted for. It was found that thermal effects significantly affect thermal load by reducing buckling load. Wave propagation in piezoelectric nanobeams, influenced by surface and thermal effects were studied in [22]. The authors pointed out the profound influence of the scale coefficient and temperature on the frequency of aforementioned nanobeams. An extension to temperature dependent properties is provided in [23, 24], with a special emphasis of the transient heating in the form of sinusoidal pulse and ramp-type heating on the vibrations of the nanobeam. It was found that the pulse width significantly affects mechanical response of the 2
ACCEPTED MANUSCRIPT
nanobeam. Other contributions to the vibration analysis of FG nanobeams can be found in [25]. Effect of electrostatic force and thermal shock on capacitive FG nanocantilevers was studied in [26]. Finally, it should be pointed out that thermal properties of FG materials are highly affected by the gradient pattern, resulting in a variety of temperature distributions, [27, 28], thus supporting the need for the reliable non-local thermomechanical model. A standard procedure to predict a material property of a nanocomposite, is based on the classical rule of mixtures. For some property P , there is: P F G = P a V a + P b Vb ,
where Va + Vb = 1.
(1)
Above, indexes a, b denote components a and b, respectively, while V represent the volume fraction. The volume fraction of the component is then specified to vary by some functional dependence, for example linear [5, 29], power law [30, 6, 3], Flory-Rehner equation [31]. Sometimes, CNT properties are scaled in order to better match experimentally noticed behaviour, [5, 19, 29]. The scaling parameter accounts for various defects in the reinforcing material or non-perfect load transfer at the CNT-matrix interface. The same approach is pursued in the present paper as well, although the paper does not set any restrictions on the particular form of functional dependence of a property. The novel aspects of this paper can be stated as follows. Starting from the strict continuum mechanics framework, augmented by the Eringen nonlocal theory [32], nonlocal thermoelastic Timoshenko beam governing equations for the bending moment, axial and shear forces are developed. However, to accommodate different nonlocal effects due to normal and shear stresses as usually assumed in Timoshenko nanobeams [33], an extension of the original formulation proposed by Eringen to the anisotropic form is now proposed. In contrast to other papers in the field, solution is offered for the case of arbitrary temperature field; usual approach assumes a homogeneous temperature field in the whole beam. Such nonisothermal and nonlocal Timoshenko formulation suitable for FG nanocomposite materials with arbitrary temperature fields was not previously addressed in the literature. Additionally, usual approach relies on the formulation of the set of ordinary differential equation (ODE) for the axial and transverse displacements as well as for the rotation angle. Unfortunately, such approach in the problems at hand leads to highly coupled set of ODE and represent an unnecessary complication. We show that much simpler way to solve the problem is determining the 3
ACCEPTED MANUSCRIPT
forces and moments first and subsequent solving for the unknown displacements and remaining integration constants. Such a simplification is most welcomed since benchmarks introduced by means of examples shows that resulting displacement fields provide rather lengthy form even for the simplest cases. The volume of the paper is divided into three main sections. In the second section, a rigorous nonlocal continuum framework is set-up. Usage of such results is performed in the third section to obtain Timoshenko nanobeam governing equations. Example section provides new benchmarks and illustrate main results of the paper. Finally, the conclusion section completes the paper and points out the most significant results. 2. Nonlocal continuum mechanics 2.1. Nonlocal stress tensor In order to introduce small-size effects, we pursue the approach introduced by Eringen [34]. Let us consider the current configuration of a body B ⊂ R at time t. It is assumed that in each point x of the body, the balance of momentum is valid: ¯ + f = ρu, ¨ divσ (2) ¯ is the nonlocal stress tensor at the point x, ρ(x, t) is the mass where σ density, u(x, t) is the displacement vector and f (x, t) are volume forces. According to [32], the nonlocal stress tensor has the form: ¯ σ(x) =
Z
α (|x0 − x| , τNL ) σ(x0 )dV (x0 ).
(3)
B
Above, α (|x0 − x| , τNL ) is the nonlocal modulus and σ(x0 ) is the classical stress tensor at point x0 . Usually, τNL = c/l, the nonlocal parameter c = e0 a, where e0 is a materially dependent constant, a is the internal characteristic length like granular distance or lattice parameter and l is the external characteristic length like crack length or similar. In that way, the stress state at x depends on stress states in all other points x0 ∈ B. Regarding the nonlocal modulus α, it can be shown that α is a Green’s function of a linear differential operator L, i.e. [32]: Lα (|x0 − x| , τNL ) = δ (|x0 − x|) .
4
(4)
ACCEPTED MANUSCRIPT
The above form includes the Dirac delta function δ ensuring that the classical (local) elasticity is encompassed in the theory for the choice c = 0, see [32] for ¯ = σ. Moreover, more details. Application of this result to Eq. (3) gives Lσ in the case when L is a differential operator with constant coefficients, body forces are not significant and inertial loads can be neglected, the classical equation of equilibrium is conveniently obtained, divσ = 0. At this point some constitutive law defining the nonlocal modulus should be introduced. Again, based on the short review of such models in [32], for this paper a suitable choice is: q
˙ 0 − x) (x0 − x)(x 1 K ( ), α (|x0 − x| , c) = 0 2πc2 c
(5)
where K0 is the modified Bessel function. For such a choice, the linear differential operator L can be shown to be [32]: L = 1 − c2 ∇2 ,
(6)
¯ and the classical stress thus providing the relation between the nonlocal σ tensor σ: ¯ = σ. (1 − c2 ∇2 )σ (7) Therefore, the nonlocal stress tensor is influenced by the scalar nonlocal parameter c that leads to the isotropic nonlocal behaviour. At this point we generalize this statement, allowing anisotropic nonlocal behaviour through the new specific tensorial form: h
i
¯ = σ, 1 − c1 ∇2 σ
(8)
¯ = σ + c1 σ1 . σ
(9)
or For the specific choice c1 = c2 I, where I represents Kronecker’s tensor, the classical Eringen form, Eq. (7), is recovered. The anisotropic form is particularly suitable for the nonlocal analysis of Timoshenko nanobeams as it will be shown later. 2.2. Balance equations Now, with the previously introduced nonlocal stress tensor, it is necessary to consider the first and the second law of thermodynamics. These laws can 5
ACCEPTED MANUSCRIPT
be written as:
¯ : d + ρr + P ρe˙ + divq = σ q − div , ρη˙ ≥ ρr θ θ
(10)
where e(x, t) the internal energy, d is the rate of deformation tensor, η is the specific entropy per unit mass, θ is temperature, r(x, t) is the heat supply per unit mass. The heat flux vector across a surface with the normal n is denoted as q. The variable P represents the rate at which the energy is supplied to a point from its neighbourhood. Note that balance of energy at the global R level implies that this rate must vanish, i.e. PdV = 0. Consequently, the Ω
balance of energy can be written in the global form as: Z
(ρe˙ + divq) dV =
Ω
Z
¯ : d + ρr) dV, (σ
(11)
Ω
while the second law of thermodynamics can be rephrased as [33]: R
¯ : ε˙ − σ : ε˙ − σ1 ∗ ∇(1) ε˙ dV = 0 σ
Ω
R Ω
(12)
− qθ grad θ dV ≥ 0.
2.3. Constitutive model The material at hand is assumed to be of the thermoelastic nature, with the nonlocal influence on the stress state in the point. Therefore, the internal energy takes the form e = e(ε, ∇ε, η), where ∇ε reflects nonlocality. At this point we do not make any assumption on the functional dependence of a material property M , so mechanical properties are considered to be given in a general form as functions of three spatial coordinates x = x(x, y, z), where x ∈ B. Regarding thermal behaviour, all properties are considered to be temperature dependent, so M = M (x, θ) = M (x, y, z, θ). For such a case, the Helmholtz free energy function ψ = e − ηθ takes the form: ¯ i ∗ ∇ε − β : ε(θ − θ0 )− ρψ(ε, ∇ε, θ) = 21 ε : C : ε + 12 ∇ε ∗ C 1 − 2θ cε (θ − θ0 )2 ,
(13)
where quantities in these expressions are: C=
∂ 2ψ , ∂ε2
¯i = C
∂ 2ψ , ∂(∇i ε)2
−β =
6
∂ 2ψ , ∂ε∂θ
1 ∂ 2ψ − cε = . θ ∂θ2
(14)
ACCEPTED MANUSCRIPT
Above, cε is the specific heat at constant strain, β is the second order tensor describing anisotropic volumetric expansion due to thermal effects, C is the ¯ i is a higher order fourth order tensor describing local elastic behaviour and C tensor describing elastic effects due to strain gradients. In the case of isotropic volumetric heat expansion, volumetric expansion is reduced to: β = 3καI,
(15)
where κ is the bulk modulus and α is the coefficient of linear thermal expansion (CTE). The anisotropic tensor describing nonlocal parameters introduced in Eq. (8) is assumed to have the specific form: c2 b2 b2 2 c1 = b c2 b2 , b2 b2 c2
(16)
thus assuming the different nonlocal parameters for normal (c) and shear (b) stresses. With the above definition, the following constitutive relationship will be used in this work: ¯ i = c1 · C = c1im Cmjkl . C
(17)
Applied to the classical elasticity, this gives; h
i
¯ i = c2 Iim + b2 (1i 1m − Iim ) [λImj Ikl + µ(Imk Ilj + Iml Ijl )] , C h
(18)
iT
where λ, µ are Lame constants and 1T = 1 1 1 . 2.4. Timoshenko beam kinematics In this section, previously developed continuum mechanics framework is specialized to the nanobeams kinematics operating in the nonisothermal environment. The in-plane bending is assumed to take place in the x − y plane, where x denotes the beam’s longitudinal axis. In order to avoid unnecessary distractions, both mechanical and thermal loading are taken to be independent of the z coordinate. Therefore, in an arbitrary beam point defined by coordinates (x, y, z), the displacement field (ux , uy , uz ) is governed
7
ACCEPTED MANUSCRIPT
by the following equations: ux (x, y) = u0 (x) − ϕ(x)y,
uy (x) = v(x),
uz = 0,
(19)
where u0 (x) represents the axial displacement caused by the axial force and the thermal dilatation originating from the homogeneous part of the temperature field, v is the transverse displacement and the rotation of the beam cross-section is ϕ. Consequently, components of the strain tensor follow in the standard manner: (1)
εx (x, y) = ε = ∂x ux (x, y) = u0 (x) − ϕ(1) (x)y, γxy (x, y) = γ = ∂x uy (x) + ∂y ux (x, y) = v (1) (x) − ϕ(x),
(20)
where exponent (•) implies derivative of the respective order with respect to the longitudinal coordinate. Eq. (20)2 implies that the rotation of the crosssection ϕ(x) is not equal to the rotation of the beam longitudinal axis v (1) (x) but it differs from it for γxy (x, y). Also, we follow the standard practice in the thermomechanics of beams and additively decompose εx into the mechanical εM and the thermal part εT . With the defined Timoshenko nanobeam kinematics framework , the internal energy can be now rephrased as e = e(ε, ε(1) , γ, γ (1) , η). The Helmoltz free energy is then: ψ(ε, ε(1) , γ, γ (1) , θ) = e − ηθ. (21) With the above Legrende transformation and the balance of energy the second law of thermodynamics, Eq. (10)2 takes the form: q σ ¯ : ε˙ + τ¯ : γ˙ − ρψ˙ − ρη θ˙ − gradθ + P ≥ 0, θ
(22)
˙ ∇ε, γ, ∇γ, θ): and introducing the time rate of the Helmholtz energy ψ(ε, ˙ (¯ σ − ρ∂ε ψ) : ε˙ + (¯ τ − ρ∂γ ψ) : γ˙ − ρ (η + ∂θ ψ) θ− q (1) (1) −ρ∂ε(1) ψ ε˙ − ρ∂γ (1) ψ γ˙ − θ gradθ + P ≥ 0.
(23)
The standard arguments [35] along with Eq. (16) are now extended to account for the nonlocal behaviour and lead to the conjugacy between stress and
8
ACCEPTED MANUSCRIPT
strain: σ = ρ∂ε ψ = σM + σT c2 σ1 = ρ∂ε(1) ψ, τ = ρ∂γ ψ, b2 τ1 = ρ∂γ (1) ψ,
(24)
so the second law of thermodynamics, see Eq. (12)1 , is now: R Ω R Ω
σ ¯ ε˙ + τ¯γ˙ − σ ε˙ − τ γ˙ − c2 σ1 ε˙(1) − b2 τ1 γ˙ (1) dV = 0, (25)
− qθ grad θ dV ≥ 0,
where dV = dx da and da = dy dz. At this point, the constitutive behaviour reported in Sec. (2.3) is specialized for the presented Timoshenko beam kinematics. In the fashion of the classical rule of mixture, Eq. (1), we introduce the following form of the free energy function at x ∈ B: ψ(x, ε, ε(1) , θ) = ψa (ε, ε(1) , θ)Va (x) + ψb (ε, ε(1) , θ)Vb (x),
(26)
where ψa , ψb are Helmholtz free energy functions describing two components. Volume fractions of each component Va , Vb are position dependent. The extension to more than two components can be carried out in straightforward manner. Restricting our elaborations to in-plane bending, i.e. dropping dependence on the coordinate z, the following form of the Helmholtz free energy for each component can be used: ρψ(ε, ε(1) , θ) = 21 E(θ)ε2 + 12 c2 E(θ)(ε(1) )2 + 1 χG(θ)γ 2 + 21 b2 χG(θ)(γ (1) )2 − 2 1 cε (θ)(θ − θ0 )2 , −E(θ)α(θ)ε(θ − θ0 ) − 2θ
(27)
where E is the Young’s modulus, G is the shear modulus and χ is the shear correction factor. Alternatively, in the case when the rule of mixture is not an appropriate selection, one can use the more general form: ρψ(x, ε, ε(1) , θ) = 21 E(x, y, θ)ε2 + 12 c2 E(x, y, θ)(ε(1) )2 + 1 χG(x, y, θ)γ 2 + 21 b2 χG(x, y, θ)(γ (1) )2 − 2 1 −E(x, y, θ)α(x, y, θ)ε(θ − θ0 ) − 2θ cε (x, y, θ)(θ − θ0 )2 ,
(28)
instead of Eq. (27). The extension to out-of-plane bending leads to the somewhat more complicated form, but it is straightforward and is left as an exercise to the reader. 9
ACCEPTED MANUSCRIPT
Also, to avoid unnecessary complication in the experimental determination of non-local parameters c and b, these parameters will be treated as temperature independent. 3. Internal forces in Timoshenko beams The starting point for the analysis of internal forces in a Timoshenko beam is Eq. (25)1 . By the application of Timoshenko beam kinematics, Eq. (20), Eq. (25)1 can be transformed into: R
(1)
(1)
(¯ σ (u˙ 0 − ϕ˙ (1) y) − (σM + σT )(u˙ 0 − ϕ˙ (1) y)−
Ω
(29)
(2)
˙ −c2 σ1 (u˙ 0 − ϕ˙ (2) y) + τ¯(v˙ (1) − ϕ)− (1) 2 (2) (1) −τ (v˙ − ϕ) ˙ − b τ1 (v˙ − ϕ˙ ))dv = 0.
In the case of the rule of mixtures, Eq. (1,26,27), the stresses Eq. (24) take the following form: σM = (Va Ea + Vb Eb )ε = EFG ε σT = −(Va Ea αa + Vb Eb αb )(θ − θ0 ) = −(Eα)FG (θ − θ0 ) c2 σ1 = c2 (Va Ea + Vb Eb )ε(1) = c2 EFG ε(1) b2 τ = χ(Va Ga + Vb Gb )γ = GFG γ b2 τ1 = χb2 (Va Ga + Vb Gb )γ = χb2 GFG γ (1) ,
(30)
where effective moduli (•)FG were introduced. Alternatively, if the material properties are defined by Eqs. (27) instead of the rule of mixtures, the appropriate functional dependence must be chosen for (•)FG . Bending moments are introduced as: ¯ = −R σ M ¯ yda, A
(1)
MM = − σM yda = − EFG εyda = ϕ(1) EFG y 2 da − u0 R
R
R
A
A
A R (2)
M1 = − σ1 yda = − EFG ε(1) yda = ϕ R
R
AR
RA
A
A
A
EFG y 2 da −
R
EFG yda,
A (2) R u0 A
EFG yda,
MT = − σT yda = (Eα)FG (θ − θ0 )yda, (31)
10
ACCEPTED MANUSCRIPT
while the shear and the axial forces are: R T¯ = τ¯da,
T =
A R
τ da =
AR
T1 =
R
χGFG γda = v (1) χ GFG da − ϕχ GFG da, R
AR
A
A
A
χGFG γ (1) da = v (2) χ GFG da − ϕ(1) χ GFG da,
τ1 da =
R P¯ = σ ¯ da,
R
A
R
R
A
A
(32)
A
PM =
R
σM da =
A
P1 = PT =
R AR A
R
EFG εda = −ϕ
(1)
A
σ1 da = EFG ε(1) da = −ϕ R
R
EFG yda +
A R (2)
A
(1) u0
EFG yda +
A
R
EFG da,
A (2) R u0 A
EFG da,
σT da = − (Eα)FG (θ − θ0 )da. R
A
Note that due to the fact that the material is functionally graded, the above integrals involve position dependent moduli and CTE. Therefore, these general integrals cannot be simplified to the usual form involving the area and the first and the second moment of inertia. However, to obtain equations in the more familiar form, the following notation for stiffnesses in a cross-section are introduced: kEI =
R
EFG y 2 da kES =
R
EFG yda,
R
AR
A
A
A
kGA = χ GFG da kEA =
EFG da.
(33)
These values are not constants but rather functions of the longitudinal coordinate x. With the above notation, Eqs. (31,32) become now: (1)
MM = ϕ(1) kEI − u0 kES , (2) M1 = ϕ(2) kEI − u0 kES , T = v (1) kGA − ϕkGA = (v (1) − ϕ)kGA , T1 = v (2) kGA − ϕ(1) kGA = (v (2) − ϕ(1) )kGA , (1) PM = −ϕ(1) kES + u0 kEA , (2) P1 = −ϕ(2) kES + u0 kEA .
11
(34)
ACCEPTED MANUSCRIPT
Eq. (29) can be now recast into the form involving bending moments, shear and axial forces: RL
(2)
(1)
[(P¯ − PM − PT )u˙ 0 − c2 P1 u˙ 0 + 0 ¯ − MM − MT + b2 T1 )ϕ˙ (1) − c2 M1 ϕ˙ (2) + +(M (T¯ − T )v˙ (1) − b2 T1 v˙ (2) + (−T¯ + T )ϕ]dx ˙ = 0,
(35)
where L represents the beam length. Application of the integration by parts, for example for the transverse displacement v˙ leads to: RL 0
RL
(•)v˙ (1) dx = [(•)v] ˙ L0 − (•)(1) vdx, ˙ 0
RL
h
(2)
(•)v˙ dx = (•)v˙
(1)
iL
h
− (•) v˙
0
0
(1)
iL 0
RL
(36) (2)
+ (•) vdx. ˙ 0
The same procedure should be applied to the variables u˙ 0 and ϕ as well. In such manner, Eq. (35) takes the following form: h
(P¯ − PM − PT )u˙ 0 h
− c
2
i (1) L P1 u˙ 0 0
h
+ c
iL
2
0
RL
(1)
0 iL (1) P1 u˙ 0 0
RL 0
¯ − MM − MT + b2 T1 )ϕ˙ + (M 2
− c M1 ϕ˙
(1)
h
iL 0
+ (T¯ − T )v˙ h
− b2 T1 v˙ (1)
iL
iL 0
0
h
+ c
2
(2)
− c2 P1 u˙ 0 dx+
h
h
(1)
− (P¯ (1) − PM − PT )u˙ 0 dx−
iL
(1) M1 ϕ˙ 0
iL 0
L
R (1) (1) (1) ¯ (1) − MM ˙ − MT + b2 T1 )ϕdx− − (M
RL
0
(2)
2
− c M1 ϕdx+ ˙ 0
RL
˙ − (T¯(1) − T (1) )vdx− 0
h
(1)
+ b2 T1 v˙
iL 0
RL
(2)
RL
− b2 T1 vdx ˙ − (T¯ − T )ϕ˙ = 0. 0
0
(37) Applying the standard localization procedure, the following set of ordinary differential equations (ODEs) is obtained: u˙ 0 : v˙ : ϕ˙ :
(1) (1) (2) P¯ (1) − PM − PT + c2 P1 = 0, (2) −T¯(1) + T (1) − b2 T1 = 0, (1) (1) (2) (1) ¯ (1) − MM M + T¯ − T − MT + c2 M1 + b2 T1 = 0,
(38)
with kinematic boundary conditions for {u0 , v, ϕ} given at {0, L} or alterna12
ACCEPTED MANUSCRIPT
tively via static conditions: (1) prescribe u0 or P¯ − PM − PT + c2 P1 = 0, (1) prescribe u0 or − c2 P1 = 0, (1) prescribe v or T¯ − T + b2 T1 = 0, prescribe v (1) or − b2 T1 = 0, ¯ − MM − MT + c2 M1(1) + b2 T1 = 0, prescribe ϕ or M prescribe ϕ(1) or − c2 M1 = 0.
(39)
¯ (1) − T¯, distributed transversal Note that distributed couples are m = −M load py = −T¯(1) and distributed longitudinal load px = −P¯ (1) , see [33] for a detailed elaboration. The transversal and the longitudinal forces in a crosssection correspond to T¯ and P¯ , respectively, while the couple in a cross¯. section corresponds to M Solution of the above ODE system can be now performed in two ways. First one is a standard one and is most frequently used, see [33] for an example. The second one, as proposed in this paper leads to a more simpler system and will be favoured in the example section. Standard approach: Consider just introduced external loads, Eqs. (31,32) and Eqs. (34). The governing ODEs Eqs. (38) can be transformed by straightforward manipulation into the final form: u˙ 0 : v˙ : ϕ˙ :
(4)
(2)
(1)
−c2 ϕ(4) kES + ϕ(2) kES + c2 u0 kEA − u0 kEA − PT − px = 0, kGA b2 (v (4) − ϕ(3) ) − kGA (v (2) − ϕ(1) ) − py = 0, (4) c2 (ϕ(4) kEI − u0 kES ) + b2 kGA (−ϕ(1) + v (2) ) + kGA (ϕ − v (1) )− (2) (1) −ϕ(2) kEI + u0 kES − MT − m = 0, (40)
13
ACCEPTED MANUSCRIPT
while the accompanying boundary conditions are: (1) (3) prescribe u0 or − c2 (ϕ(3) kES − u0 kEA ) + ϕ(1) kES − u0 kEA − PT + P¯ = 0, (1) (2) prescribe u0 or ϕ(2) kES − u0 kEA = 0, prescribe v or − b2 kGA (ϕ(2) − v (3) ) + kGA ϕ − kGA v (1) + T¯ = 0, prescribe v (1) or b2 kGA (ϕ(1) − v (2) ) = 0, (3) prescribe ϕ or c2 (ϕ(3) kEI − u0 kES ) − b2 kGA (ϕ(1) − v (2) )− (1) ¯ = 0, −ϕ(1) kEI + u0 kES − MT + M (2) (1) 2 (2) prescribe ϕ or c (−ϕ kEI + u0 kES ) = 0. (41) It can be easily seen, the resulting system of ODE is significantly coupled. Alternative approach: Alternatively to the usual practice, we do not transform the problem to the form involving displacements and the angle as unknowns. Instead, we rather choose to work with forces and moments as unknowns. First, note that from Eqs. (34) follows: (1)
M1 = MM ,
T1 = T (1) ,
(1)
P1 = PM .
(42)
Applying the above result, the Eqs. (38) take the following form u˙ 0 : v˙ : ϕ˙ :
(3) (1) (1) P¯ (1) − PM − PT + c2 PM = 0, −T¯(1) + T (1) − b2 T (3) = 0, (3) (1) (1) ¯ (1) − MM + T¯ − T − MT + c2 MM + b2 T (2) = 0, M
(43)
with boundary conditions : (2) prescribe u0 or P¯ − PM − PT + c2 PM = 0, (1) (1) prescribe u0 or − c2 PM = 0, prescribe v or T¯ − T + b2 T (2) = 0, prescribe v (1) or − b2 T (1) = 0, ¯ − MM − MT + c2 M (2) + b2 T (1) = 0, prescribe ϕ or M M (1) (1) 2 prescribe ϕ or − c MM = 0.
(44)
Obviously, Eqs. (43)1,2 are decoupled and can be easily solved. Subsequently, with the shear and axial forces defined, the solution of Eq. (43)3 is straightforward.
14
ACCEPTED MANUSCRIPT
Some of unknown integration constants follow from the static boundary conditions. This is not always the case. For example, if the doubly clamped beam is analysed, not all integration constants can be determined from static boundary conditions. Unknown derivatives of displacement functions can be now easily obtained from the solution of the linear system Eqs. (34). Application of the kinematic boundary conditions is used to determine integration constants in this operation as well as the remaining integration constants (if any) from the solving of static quantities. From another point of view, the problem can be also imagined as the set of 6 ODE, Eqs. (43) and Eqs. (34)1,3,5 , six unknown functions {PM , T, MM , u0 , v, ϕ}, accompanied by the suitable set of boundary conditions Eqs. (44). 4. Examples 4.1. Functionally graded cantilever beam loaded by the non-homogeneous temperature field Let us consider a cantilever nanobeam of length L and of the rectangular cross-section HxB. The longitudinal axis is denoted as x while bending is assumed to take place in x − y plane. Young’s and shear modulus are linear functions of the coordinate y: EFG = E0 + E1 y,
GFG = G0 + G1 y,
(Eα)FG = Eα0 + Eα1 y,
(45)
where E0 , E1 , G0 , G1 are some material constants. The product (Eα)FG also assumes a linear variation along the coordinate y defined by additional material parameters Eα0 , Eα1 . The temperature field is non-homogeneous and is defined by the quadratic distribution θ = θ(y) along the height of the nanobeam: ∆θ(y) = t0 + t1 y + t2 y 2 . (46) Additionally, a case θ = θ(x, y) where the above distribution is scaled by the longitudinal coordinate x is considered as well: ∆θ(x, y) = (t0 + t1 y + t2 y 2 )x.
(47)
The beam is not mechanically loaded in any way and is clamped at the right ¯ (0) = 0 and the Eqs. (43) take the end. Consequently,P¯ (0) = T¯(0) = 0, M
15
ACCEPTED MANUSCRIPT
following form (1)
(1)
(3)
−PM − PT + c2 PM = 0, T (1) − b2 T (3) = 0, (1) (1) (3) −MM + T¯ − T − MT + c2 MM + b2 T (2) = 0.
(48)
Boundary conditions, Eqs. (39), are: at x = 0 (2) −PM − PT + c2 PM = 0, −T + b2 T (2) = 0, (2) −MM − MT + c2 MM + b2 T (2) = 0, (1) (1) −c2 PM = 0, −b2 T (1) = 0, −c2 MM = 0 at x = L u0 = 0, v = 0, ϕ = 0, (1) (1) −c2 PM = 0, −b2 T (1) = 0, −c2 MM = 0.
(49)
Case θ = θ(y) : In the case when temperature is defined by the Eq. (46), solutions of Eq. (48) are: PM = cex/c c1 − ce−x/c c2 + c3 , T = c4 + 21 be−x/b (−1 + e2x/b )c5 + 12 b2 e−x/b (−1 + ex/b )2 c6 , x x MM = c7 + 12 ce− c (−1 + e2x/c )c8 + 12 c2 e− c (−1 + ex/c )2 c9 + x 2x x 2x + 12 e− c (−c + ce c − 2xex/c )c4 − 12 b2 e− c (−c + ce c − 2xex/c )c6 .
(50)
Solving for the integration constants, it is obtained: c1 = c2 = c4 = c5 = c6 = c8 = c9 = 0, 1 c3 = BHEα0 t0 + 12 BH 3 (Eα1 t1 + Eα0 t2 ), 1 1 3 c7 = − 12 BH (Eα0 t1 + Eα1 t0 ) − 80 BH 5 Eα1 t2 ,
(51)
thus collapsing the solution to: 1 PM = c3 = BHEα0 t0 + 12 BH 3 (Eα1 t1 + Eα0 t2 ), T = 0, 1 1 MM = c7 = − 12 BH 3 (Eα0 t1 + Eα1 t0 ) − 80 BH 5 Eα1 t2 ,
(52)
i.e, both PM and MM are independent of the coordinates x while the shear 16
ACCEPTED MANUSCRIPT
force vanishes. Now, from Eqs. (34) it is: (1)
PM = −ϕ(1) kES + u0 kEA , 0 = (v (1) − ϕ)kGA , (1) MM = ϕ(1) kEI − u0 kES ,
(53)
along with boundary conditions u0 (L) = 0, v(L) = 0, ϕ(L) = 0. Integrating for the unknown functions, it is obtained: M E0 +MM E1 )x + c10 , u0 = − 12(P BH(E 2 H 2 −12E 2 ) 1
0
2
2
M E0 +PM E1 H )x v = − (6(12M + c11 + xc12 , BH 3 (−12E 2 +E 2 H 2 )) 0
ϕ=
1
2 M E0 +PM E1 H )x − (12(12M BH 3 (−12E02 +E12 H 2 ))
(54)
+ c12 ,
with integration constants: c10 =
2 2 2 (12(PM E0 +MM E1 )L) M E0 L +6PM E1 H L , c11 = 72M12BE , 2 H 3 −BE 2 h5 BH(−12E02 +E12 H 2 ) 0 1 2 M E0 +PM E1 H )L c12 = 12(12M . BH 3 (−12E02 +E12 H 2 )
(55)
Thus, the axial displacement and the bending angle are linearly distributed, while the transverse displacement takes the parabolic form. Case θ = θ(x, y) : In the case when temperature is defined by the Eq. (47), solutions of Eqs. (48) are: 1 B(12Eα0 Ht0 + Eα1 H 3 t1 + Eα0 H 3 t2 )x + cex/c c1 − ce−x/c c2 + c3 , PM = 12 T = c7 + 21 be−x/b (−1 + e2x/b )c8 + 21 b2 e−x/b (−1 + ex/b )2 c9 , 1 MM = 960 cBH 3 (−1 + ex/c )2 (1 − e−2x/c )(20Eα1 t0 + 20Eα0 t1 + 3Eα1 H 2 t2 )− 1 cBH 3 (1 + e−2x/c )(−1 + e2x/c )(20Eα1 t0 + 20Eα0 t1 + 3Eα1 H 2 t2 )+ 960 1 BH 3 (20Eα1 t0 + 20Eα0 t1 + 3Eα1 H 2 t2 )(−ce−x/c + cex/c − 2x) + c4 + 480 1 −x/c ce (−1 + e2x/c )c5 + 21 c2 e−x/c (−1 + ex/c )2 c6 + 21 e−x/c (−c + ce2x/c − 2ex/c x)c7 − 2 1 2 −x/c be (−c + ce2x/c − 2ex/c x)c9 , 2 (56)
17
ACCEPTED MANUSCRIPT
Solving for the integration constants, it is obtained: 2
α1 H t1 +Eα0 H c1 = − BH(12Eα0 t0+E 12(1+eL/c )
c2 c3 c4 c5 c6 c7
2t ) 2
,
BH(12Eα0 t0 +Eα1 H 2 t1 +Eα0 H 2 t2 ) 1 B(12Eα0 Ht0 + Eα1 H 3 t1 − 12 12(1+eL/c ) 1 (12BEα0 ht0 x + BEα1 H 3 t1 x + BEα0 H 3 t2 x), 12 3 2 L/c )+x+eL/c x)) α1 H t2 )(−c+ce − BH (20Eα1 t0 +20Eα0 t1 +3E , 240(1+eL/c ) 1 3 2 BH (20Eα1 t0 + 20Eα0 t1 + 3Eα1 H t2 ), 240 BH 3 (−1+eL/c )(20Eα1 t0 +20Eα0 t1 +3Eα1 H 2 t2 ) − , 240(1+eL/c )
= = = = = = c8 = c9 = 0,
+ Eα0 H 3 t2 ),
(57) thus collapsing the solution Eqs. (48) to: −x/c
2
2
L/c
2x/c
x/c
L/c+x/c
H t2 )(−ce +ce −2e x−2xe ) PM = − BHe (12Eα0 t0 +Eα1 H t1 +Eα012(1+e L/c ) T = 0, 3 −x/c 2 t )(−ceL/c +ce2x/c −2xex/c −2xeL/c+x/c ) 2 MM = BH e (20Ea1t0 +20Ea0t1 +3Ea1H . 240(1+eL/c )
(58) Differently from the previous case, PM and MM are now dependent on the coordinate x, while the shear force vanishes once more. Again, from Eqs. (34) and with boundary conditions u0 (L) = 0, v(L) = 0, ϕ(L) = 0,it is obtained: L−x c +c2 (−ex/c )+x2 (eL/c +1))(20E0 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 H 2 t1 )−E1 H 2 (20Eα0 t1 +3Eα1 H 2 t2 +20Ea1t0 )) 20(eL/c +1)(12E02 −E12 H 2 ) L−x 2 3 L/c 2 2 x/c c (−c (c+x)e −c e (x−c)+ 3 x (e +1))5E1 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 H 2 t1 )−3E0 (20Eα0 t1 +3Eα1 H 2 t2 +20Eα1 t0 ) 5(eL/c +1)(12E02 −E12 H 2 ) L−x 2 2 x/c 2 L/c 2 x(−c e c −c e +x (e +1))(5E1 (Eα0 H t2 +12Eα0 t0 +Eα1 H 2 t1 )−3E0 (20Eα0 t1 +3Eα1 H 2 t2 +20Eα t0 )) 12 11 5(eL/c +1)(12E02 −E12 H 2 ) L−x 2 2 x/c 2 L/c 2 2 2 c (−c e +c (−e )+x (e +1))(5E1 (Eα0 H t2 +12Eα0 t0 +Eα1 H t1 )−3E0 (20Eα0 t1 +3Eα1 H t2 +20Eα1 t0 )) 12 5(eL/c +1)(12E02 −E12 H 2 )
u0 =
(−c2 e
+ c10 ,
v=
+
+c x+c ,
ϕ=
+c ,
(59) with integration constants: c10 =
(c−L)(c+L)(20E0 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 H 2 t1 )−E1 H 2 (20Eα0 t1 +3Eα1 1H 2 t2 +20Eα1 t0 )) , 240E02 −20E12 H 2 −L/2c
L/2c
c11 =
(3c3 eL/2c −e−L/2c −3c2 L+2L3 )(5E1 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 H 2 t1 )−3E0 (20Eα0 t1 +3Eα1 H 2 t2 +20Eα1 t0 )) e
c12 =
+e
180E02 −15E12 H 2 (c−L)(c+L)(5E1 (Eα0 H 2 t2 +12Eα0 t0 +Eα1 a1H 2 t1 )−3E0 (20Eα0 t1 +3Eα1 H 2 t2 +20Eα1 t0 )) 60E02 −5E12 H 2
(60) Therefore, if these two cases with different temperature distributions are
18
.
,
ACCEPTED MANUSCRIPT
compared, it can be said that in the beam that is not loaded by mechanical forces, shear forces do not occur. As expected, longitudinal temperature variation leads toward more complex solutions. In particular, the axial displacement in the case of θ = θ(y) is distributed linearly, while in the second case significantly more complex distribution involving exponential functions take place. 4.2. Simply supported functionally graded beam loaded by non-homogeneous temperature field and mechanical quantities Let us consider a simply supported nanobeam of L = 100 nm and of the rectangular cross-section H = 20 nm, B = 10 nm. As in the previous example, the beam longitudinal axis is denoted with x and deformation is assumed to take place in x − y plane. At the left end, the beam elongation is constrained in all directions, ie. u(0) = v(0) = 0, while at the right end beam dilatation along x axis is free, so v(L) = 0. Shear correction factor is χ = 5/6. The nanobeam is made of Poly(methyl methacrylate) (PMMA) matrix. Young’s modulus of PMMA is EPMMA = 2.5 GPa, shear modulus GPMMA = 0.93 GPa, specific density ρPMMA = 1190 kg/m3 and the coefficient of thermal expansion αPMMA = 70 · 10−6 K−1 . As a reinforcement, single walled carbon nanotubes (CNT) are used. CNTs Young’s modulus is ECNT = 1 TPa, GCNT = 0.5 TPa and αCNT = −1.2 · 10−5 K−1 . Note that the negative sign of CTE describes contraction of the CNT with the temperature increase. Nanotubes are dispersed with approximately linear variation along CNT the beam height, as depicted in Fig. 1, i.e. 30 CNTs in total. Therefore, CNT volume fraction is: PMMA
CNT PMMA
Figure 1: The cross section of the simply supported nanobeam
VCNT = 2
|y| max , V H CNT
19
(61)
ACCEPTED MANUSCRIPT
max where VCNT represent maximal volume fraction of CNTs at the y = ±H/2. max In this case VCNT = 0.1. Consequently, according to the Eq. (1), PMMA volume fraction must be:
VPMMA = 2
H − |y| Vmax . H
(62)
Three loading cases are considered: Case a. The beam is thermally loaded by the homogeneous temperature field, ∆θ = aa = 30 K. Case b. The beam is thermally loaded by the linear variation of temperature field along the longitudinal axis, ∆θ = ab x = 0.3x K. Case c. The beam is thermally loaded by the linear variation of temperature field along the longitudinal axis, ∆θ = ac = 0.3x K and by the downward distributed load q = 5 nN/nm. According to the above defined loading, it is px = 0, m = 0. In Cases a and b, py = 0, while in the Case c it is py = −q. Eventually, Eqs. (43) collapse to: u˙ 0 : v˙ : ϕ˙ :
(1) (3) P¯ (1) − PM + c2 PM = 0, −py + T (1) − b2 T (3) = 0, (1) (3) −MM + T¯ − T + c2 MM + b2 T (2) = 0.
(63)
¯ (0) = M ¯ (L) = 0 and due Since the beam is simply supported, P¯ (L) = 0, M to boundary conditions, Eqs. (39), it is : at x = 0 u0 = 0, v = 0 (2) −MM − MT + c2 MM + b2 T (2) = 0, (1) (1) −c2 PM = 0, −b2 T (1) = 0, −c2 MM = 0, at x = L (2) −PM − PT + c2 PM = 0, v = 0, (1) −MM − MT + c2 MM + b2 T (1) = 0, (1) (1) −c2 PM = 0, −b2 T (1) = 0, −c2 MM = 0. 20
(64)
ACCEPTED MANUSCRIPT
Upon determination of PM , T, MM , the displacements and the rotation can be evaluated. Temperature distribution does not depend on the longitudinal coordinate, so ∂x kθA = 0, ∂x kθS = 0. Consequently, system of equations provided by Eqs. (34) collapses to: (1)
MM = ϕ(1) kEI − u0 kES , T = v (1) kGA − ϕkGA = (v (1) − ϕ)kGA , (1) PM = −ϕ(1) kES + u0 kEA .
(65)
With the above equations at hand, for each of the three described loading cases, the solutions are obtained as follows. 4.2.1. Case a. The thermal axial force and the thermal bending moment, according to the Eqs. (31)4 ,(32)7 , are PT =
R
σT da = − (Eα)FG (θ − θ0 )da = 45.6 nN, R
A
A
MT = − σT yda = (Eα)FG (θ − θ0 )yda = 8.1 nNnm. R
R
A
A
(66)
With the above values, solution of the ODE system Eqs. (70) is: PM = −PT = −45.6 nN, T = 0, MM = −MT = −8.1 nNnm.
(67)
Consequently, both PM and MM do not change along the beam axis, while shear forces do not occur. Now, displacements and rotation fields can be evaluated from Eqs. (65): ES MT x = −23.98 · 10−5 x, u = − KKEIEAPTK+K 2 EI −K ES
v= ϕ=
(KEA LMT +KES LPT )x−(KEA MT +KES PT )x2 = 1.14 · 10−6 x 2 2(KEA KEI −KES ) (L−2x)(KEA MT +KEs PT ) = 1.14 · 10−6 − 2.272 · 10−8 x. 2 2(KEA KEI −KES )
− 1.136 · 10−8 x2 ,
(68) These results are presented in Figs. 2-4. Evidently, nonlocal effects play no role in the case of homogeneous temperature field. Please note that BernoulliEuler and Timoshenko formulations give identical results for a such case. Negative sign in the axial displacement field indicate that the effect of thermal 21
ACCEPTED MANUSCRIPT
contraction of carbon nanotubes dominates - leading toward contraction of the nanobeam. 0.000 -0.005 uHnmL
-0.010 -0.015 -0.020 -0.025
0
20
40
60
L HnmL
80
100
Figure 2: Longitudinal displacements of the simply supported nanobeam in the case of the constant temperature field
0.00003 0.000025 0.00002 vHnmL 0.000015 0.00001 5. ´ 10 -6 0
0
20
40
60
L HnmL
80
100
Figure 3: Transversal displacements of the simply supported nanobeam in the case of the constant temperature field
4.2.2. Case b. In this case, the thermal axial force and the thermal bending moment Eqs. (31)4 ,(32)7 , are not constant along the nanobeam: PT =
R
σT da = − (Eα)FG (θ − θ0 )da = PT0 x, R
A
A
MT = − σT yda = (Eα)FG (θ − θ0 )yda = MT0 x, R
R
A
A
22
(69)
ACCEPTED MANUSCRIPT
1.5 ´ 10 -6 1. ´ 10 -6 5. ´ 10 -7
j HradL
0 -5. ´ 10 -7 -1. ´ 10 -6 -1.5 ´ 10 -6
0
20
40
60
80
L HnmL
100
Figure 4: Rotation angle of the simply supported nanobeam in the case of the constant temperature field
where
max max b PT0 = − BHa + 2E2 α2 ), − E2 α2 VCNT (E1 α1 VCNT 2 2 max BH VCNT ab MT0 = (E1 α1 − E2 α2 ). 12
(70)
Consequently, both PT and MT along the beam axis vary in the linear manner. With the above values, solution of the ODE system Eqs. (70) is:
L
L
x
x
PT0 −xe c −ce c − c +ce c −x
PM = MM =
MT0
L e c +1 L L x x −xe c −ce c − c +ce c −x
(71) ,
L
e c +1
T = 0. Once more, shear forces do not occur. Having in mind that KEA , KES and KEI are not functions of the longitudinal coordinate x, displacements and rotation fields can be evaluated from Eqs. (71):
u=
L
L
x
x
L
2 e c +1
h
v=
L
2c2 e c −2c2 e c − c +2c2 −2c2 e c +x2 e c +x2 (kEI PT0 +kES MT0 )
L
x
2 (kEA kEI −kES )
L
L
L
6L e c +1
ϕ=
i
L
6c3 (Le c −L+Le c −2xe c −LeL/c−x/c +2x)+(L3 xe c +L3 x−Lx3 e c −Lx3 ) (kEA MT0 +kES PT0 )
L
x
L
x
2 (kEA kEI −kES )
L
L
−12c3 e c +12c3 +6c2 Le c +6c2 Le c − c +L3 e c +L3 −3Lx2 e c −3Lx2 (kEA PM0 +kES PT0 )
L
6L e c +1
2 (kEA kEI −kES )
. (72)
23
ACCEPTED MANUSCRIPT
0 -10 Pm HnNL
-20 -30
c=1
-40
c=20
-50
c=1000
-60 0
20
40
60
80
100
LHnmL
Figure 5: Variation of Pm along the simply supported nanobeam under the linear temperature field for increasing values of nonlocal parameter c 0
Mm HnNnmL
-2 -4 c=1
-6
c=20
-8
c=1000
-10 -12
0
20
40
60
80
100
LHnmL
Figure 6: Variation of Mm along the simply supported nanobeam under the linear temperature field for increasing values of nonlocal parameter c 0.000
u HnmL
-0.005 c=1
-0.010
c=20 c=1000
-0.015 0
20
40
60
80
100
LHnmL
Figure 7: Longitudinal displacements of the simply supported nanobeam under the linear temperature field for increasing values of the nonlocal parameter c
24
ACCEPTED MANUSCRIPT
0.00002
vHnmL
0.000015 c=1
0.00001
c=20
5. ´ 10-6
c=1000
0
0
20
40
60
80
100
LHnmL
Figure 8: Transversal displacements of the simply supported nanobeam under the linear temperature field for increasing values of the nonlocal parameter c
j HradL
5. ´ 10-7 0
c=1 c=20
-5. ´ 10-7
c=1000
-1. ´ 10-6 0
20
40
60
80
100
LHnmL
Figure 9: Rotation angle of the simply supported nanobeam under the linear temperature field for increasing values of the nonlocal parameter c
25
ACCEPTED MANUSCRIPT
Note that PM and MM are described by the same class of curves and differ by a factor PT0 /MT0 , compare Figs. 5,6. Nonlocal influence on the transverse displacement v is manifested through a slight increase of the maximal displacement and shift of the extremum away from from the midpoint span. 4.2.3. Case c. Following the same procedure as before, the solution can be readily obtained. The thermal axial force and the thermal bending moment take same values as in the Case b, PT = PT0 x, MT = MT0 x. Now, the solution of the ODE system, Eqs. (70) is:
PM =
x
L
L
x
PT0 −xe c −ce c − c +ce c −x
T =− MM =
py
,
L e c +1 L x L L L−x L2 e b +L2 +2bLe b −2Lx(e b +1)−4be b −2bLe b b +4b
e 2L(
L/b +1
2L
1
L e b +1
L − xc c
2L e c −1
L
L
,
) L
x
L
L
L
2L
L
py e c L2 c(e b + e c + e b + c + 1)
L
L
L
2L
+e L2 c(e b + e b + c + e c + 1) + 2Lc2 (e b − e c − e b + c + 1) L L L L x L L L x L L +4bc(e c − c (e b − e b + c + e c − 1) − e c (e b + e c − e b + c − 1)) L
+2b(2x − L)(e b + e
L
2L c
L
2L c
−eb+
L
L
x
L
− 1) + xL(x − L)(e b − e L
x
L
x
2L c
L
−eb+
2L c
+ 1)
x
+2MT0 L c(1 − e c )(e c − c + e b + c − c − e b + c − e c ) L
+x(e b − e
2L c
L
−eb+
2L c
+ 1)
. (73)
Distribution of PM is identical to the result obtained in Case b. Unfortunately, resulting functions that describe displacements and the rotation take rather lengthy form and due to economy in space, the formal results are not presented. Instead, these are presented in the graphical form, thus illustrating main points. As visible from Fig. 10 and Eqs. (73)2 , distribution of T is insensitive of the nonlocal parameter c. From Fig. 11 it can be noticed that the influence of the nonlocal parameter b is much more pronounced than the nonlocal parameter c. The bending moment MM is pretty much insensitive of the nonocal parameter b. The axial displacement at the midpoint of the nanobeam, Fig. 12, does not depend on b too, while transversal displacements and rota-
26
ACCEPTED MANUSCRIPT
tion angle, Figs. 13 and 14, depend on both nonlocal parameters. Sensitivity with respect to both b and c is approximately equal when the transversal displacement is considered, while influence on nonlocal parameters on the rotation angle manifest different trends. 150
b
100
50
0 4 2 T 0
50 100
0
c
150
Figure 10: T at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
0
50 c 100 150 6000 5500 5000 M M 4500 0 50
100
b
150
Figure 11: MM at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
5. Conclusions In the paper at hand a model for nanoscale thermomechanical analysis of beams was presented. In order to properly describe effects occurring at small scales, the extended version of Erignen’s original nonlocal formulation was proposed. The nanobeam model was based on Timoshenko formulation. 27
ACCEPTED MANUSCRIPT
- 0.011 - 0.012 - 0.013 - 0.014 u
0 50 c
0 50
100
100 b 150 150 0
50 c
100 150
150
100 b
50
0
- 0.011 - 0.012 - 0.013 - 0.014 u
Figure 12: The axial displacement at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
28
ACCEPTED MANUSCRIPT
- 0.9 v 1.0 - 1.1 150
150 100
100 c
50
50
b
00
- 0.9 - 1.0
150
- 1.1 150
100 c
v
100
50
50
b
0 0
Figure 13: The transversal displacement at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
29
ACCEPTED MANUSCRIPT
c
0
50
100 150 0.0002 j 0.0001 0.0000 150
100
50
b c 150
100
50
0
0
0.0002 0.0001 0.0000 j 150
100 b
50
0
Figure 14: The rotation angle at the midpoint of the simply supported nanobeam for various values of nonlocal parameters b, c - Case c
30
ACCEPTED MANUSCRIPT
The approach was verified on two examples, thus establishing two new benchmarks in the field. The main findings point out necessity of employment of the Timoshenko formulation instead of the Bernoulli-Euler formulation if the loading, either mechanical or thermal, takes other form than the most simplest. Thus, although the Bernoulli-Euler formulation or isothermal models are simpler to implement, at the same time they offer reduced possibilities of mechanical structural modelling at the nano level. In that way, the present formulation surpasses previous findings providing more general solutions and wider field of application. Moreover, the strict thermodynamical setting assures that second law of thermodynamics is fulfilled for all possible scenarios, hence providing the physically sound formulation. 6. Acknowledgements This work has been fully supported by Croatian Science Foundation under the project no. 6876 - Assessment of structural behaviour in limit state operating conditions. This support is gratefully acknowledged. References References [1] Liew K, Lei Z, Zhang L. Mechanical analysis of functionally graded carbon nanotube reinforced composites: A review. Composite Structures 2015;120(0):90 –7. [2] Kiani K. Longitudinal and transverse instabilities of moving nanoscale beam-like structures made of functionally graded materials. Compos Struct 2014;107:610–9. [3] Nazemnezhad R, Hosseini-Hashemi S. Nonlocal nonlinear free vibration of functionally graded nanobeams. Compos Struct 2014;110(1):192–9. [4] Hosseini-Hashemi S, Nazemnezhad R, Bedroud M. Surface effects on nonlinear free vibration of functionally graded nanobeams using nonlocal elasticity. Appl Math Modell 2014;38(14):3538–53. [5] Ke LL, Yang J, Kitipornchai S. Nonlinear free vibration of functionally graded carbon nanotube-reinforced composite beams. Compos Struct 2010;92(3):676–83. 31
ACCEPTED MANUSCRIPT
[6] Simsek M. Non-linear vibration analysis of a functionally graded Timoshenko beam under action of a moving harmonic load. Compos Struct 2010;92(10):2532–46. [7] Yas M, Heshmati M. Dynamic analysis of functionally graded nanocomposite beams reinforced by randomly oriented carbon nanotube under the action of moving load. Appl Math Modell 2012;36(4):1371–94. [8] Simsek M, Yurtcu H. Analytical solutions for bending and buckling of functionally graded nanobeams based on the nonlocal Timoshenko beam theory. Compos Struct 2013;97:378–86. [9] Eltaher M, Khairy A, Sadoun A, Omar FA. Static and buckling analysis of functionally graded Timoshenko nanobeams. Appl Math Comput 2014;229:283–95. [10] Eltaher M, Emam S, Mahmoud F. Static and stability analysis of nonlocal functionally graded nanobeams. Compos Struct 2013;96:82–8. [11] Barretta R, Feo L, Luciano R, Marotti de Sciarra F. Variational formulations for functionally graded nonlocal Bernoulli-Euler nanobeams. Compos Struct 2015;129:80–9. [12] Marotti de Sciarra F, Canadija M, Barretta R. A gradient model for torsion of nanobeams. C R Mec 2015;343(4):289–300. [13] Barretta R, Feo L, Luciano R. Torsion of functionally graded nonlocal viscoelastic circular nanobeams. Composites Part B: Engineering 2015;72:217–22. [14] Erisken C, Kalyon D, Wang H. Functionally graded electrospun polycaprolactone and β -tricalcium phosphate nanocomposites for tissue engineering applications. Biomaterials 2008;29(30):4065–73. [15] Filiz S, Aydogdu M. Wave propagation analysis of embedded (coupled) functionally graded nanotubes conveying fluid. Compos Struct 2015;132:1260–73. [16] Friedrich K, Zhang Z, Schlarb A. Effects of various fillers on the sliding wear of polymer composites. Compos Sci Technol 2005;65(15-16 SPEC. ISS.):2329–43. 32
ACCEPTED MANUSCRIPT
[17] Kwon H, Bradbury C, Leparoux M. Fabrication of functionally graded carbon nanotube-reinforced aluminum matrix composite. Adv Eng Mater 2011;13(4):325–9. [18] Canadija M, Barretta R, de Sciarra M. A gradient elasticity model of Bernoulli-Euler nanobeams in non-isothermal environments. submitted ????;. [19] Shen HS. Postbuckling of nanotube-reinforced composite cylindrical shells in thermal environments, part I: Axially-loaded shells. Compos Struct 2011;93(8):2096–108. [20] Shen HS, Xiang Y. Nonlinear vibration of nanotube-reinforced composite cylindrical shells in thermal environments. Comput Meth Appl Mech Eng 2012;213-216:196–205. [21] Khater M, Eltaher M, Abdel-Rahman E, Yavuz M. Surface and thermal load effects on the buckling of curved nanowires. Engineering Science and Technology, an International Journal 2014;17(4):279 –83. [22] Zhang YW, Chen J, Zeng W, Teng YY, Fang B, Zang J. Surface and thermal effects of the flexural wave propagation of piezoelectric functionally graded nanobeam using nonlocal elasticity. Comput Mater Sci 2015;97:222–6. [23] Zenkour A, Abouelregal A. Nonlocal thermoelastic nanobeam subjected to a sinusoidal pulse heating and temperature-dependent physical properties. Microsyst Technol 2014;:1–10. [24] Zenkour A, Abouelregal A, Alnefaie K, Abu-Hamdeh N, Aifantis E. A refined nonlocal thermoelasticity theory for the vibration of nanobeams induced by ramp-type heating. Appl Math Comput 2014;248:169–83. [25] Alshorbagy A. Temperature effects on the vibration characteristics of a functionally graded thick beam. Ain Shams Engineering Journal 2013;4(3):455–64. [26] JafarSadeghi-Pournaki I, Zamanzadeh M, Madinei H, Rezazadeh G. Static pull-in analysis of capacitive FGM nanocantilevers subjected to thermal moment using Eringen’s nonlocal elasticity. Int J Eng Trans A 2014;27(4):633–42. 33
ACCEPTED MANUSCRIPT
[27] Bafekrpour E, Simon G, Habsuda J, Naebe M, Yang C, Fox B. Fabrication and characterization of functionally graded synthetic graphite/phenolic nanocomposites. Mater Sci Eng, A 2012;545:123–31. [28] Bafekrpour E, Simon G, Yang C, Chipara M, Habsuda J, Fox B. Functionally graded carbon nanofiber-phenolic nanocomposites for sudden temperature change applications. Polymer (United Kingdom) 2013;54(15):3940–8. [29] Mehrabadi S, Aragh B, Khoshkhahesh V, Taherpour A. Mechanical buckling of nanocomposite rectangular plate reinforced by aligned and straight single-walled carbon nanotubes. Composites Part B: Engineering 2012;43(4):2031–40. [30] Shen HS, Noda N. Postbuckling of pressure-loaded fgm hybrid cylindrical shells in thermal environments. Compos Struct 2007;77(4):546–60. [31] Ahankari S, Kar K. Processing of styrene butadiene rubber-carbon black nanocomposites with gradation of crosslink density: Static and dynamic mechanical characterization. Mater Sci Eng, A 2008;491(1-2):454–60. [32] Eringen A. On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves. J Appl Phys 1983;54(9):4703–10. [33] Marotti De Sciarra F, Barretta R. A gradient model for Timoshenko nanobeams. Physica E 2014;62:1–9. [34] Eringen A. Nonlocal Continuum Field Theories. Springer Verlag; 2002. [35] Coleman B, Gurtin M. Thermodynamics with internal state variables. J Chem Phys 1967;47:597–613.
34