Second strain gradient theory in orthogonal curvilinear coordinates: Prediction of the relaxation of a solid nanosphere and embedded spherical nanocavity

Applied Mathematical Modelling 76 (2019) 669–698

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Second strain gradient theory in orthogonal curvilinear coordinates: Prediction of the relaxation of a solid nanosphere and embedded spherical nanocavity Farzaneh Ojaghnezhad a,∗, Hossein M. Shodja b,c a

School of Engineering, Alzahra University, Vanak St., 1993891176, Tehran, Iran Department of Civil Engineering, Sharif University of Technology, P.O. Box 11155-4313, Tehran, Iran c Institute for Nanoscience and Nanotechnology, Sharif University of Technology, P.O. Box 11155-9161, Tehran, Iran b

a r t i c l e

i n f o

Article history: Received 27 September 2018 Revised 10 June 2019 Accepted 17 June 2019 Available online 20 June 2019 Keywords: Second strain gradient theory Orthogonal curvilinear coordinates Surface effect Nanosphere Nanocavity Relaxation

a b s t r a c t In this paper, Mindlin’s second strain gradient theory is formulated and presented in an arbitrary orthogonal curvilinear coordinate system. Equilibrium equations, generalized stressstrain constitutive relations, components of the strain tensor and their first and second gradients, and the expressions for three different types of traction boundary conditions are derived in any orthogonal curvilinear coordinate system. Subsequently, for demonstration, Mindlin’s second strain gradient theory is represented in the spherical coordinate system as a highly-practical coordinate system in nanomechanics. Second strain gradient elasticity have been developed mainly for its ability to capture the surface effects in the presence of micro-/nano- structures. As a numeric illustration of the theory, the surface relaxation of spherical domains in Mindlin’s second strain gradient theory is considered and compared with that in the framework of Gurtin–Murdoch surface elasticity. It is observed that Mindlin’s second strain gradient theory predicts much larger value for the radial displacement just near the surface in comparison to Gurtin–Murdoch surface elasticity. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Recently, there has been a flurry of interest in such unusual forms as nanowires, nanotubes, and nano-particles of commonplace materials like metals, semiconductors, insulators, and organic compounds. In view of their vast possible applications in electronics, energy conversion, optics, chemical sensing, cancer therapy, and drug delivery, among other fields, consideration of the mechanics and physics of these forms of nano-structures is inevitable. The most important key features in these structures are their nanometer scale in two/three dimensions as well as their special symmetrical shapes. It is well-known that the traditional continuum theories are inadequate in treating the mechanical aspects of nano-scale structures. Due to this inadequacy, some authors have considered lattice models with the surface relaxation taken into account to examine the size-dependent elastic constants of a nanofilm [1,2]. In the context of continuum approaches, however, resorting to augmented continuum theories seems to be remedial [2,3]. Lazar [4] deals with strain gradient elasticity of defects to give a non-singular dislocation continuum theory. Lazar and Agiasofitu [5] provide fundamental quantities in generalized elasticity and dislocation theory of crystals. The development of higher order continuum theories such as strain ∗

Corresponding author. E-mail address: [email protected] (F. Ojaghnezhad).

https://doi.org/10.1016/j.apm.2019.06.021 0307-904X/© 2019 Elsevier Inc. All rights reserved.

670

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

gradient elasticity has been brought into focus, mainly in the period of about 1960–1975. In first strain gradient theory, Toupin [6] assumed that the potential energy density function of the material depends on both the second order strain tensor and its first gradient. The correspondence between first strain gradient theory and the atomic structure of the material is exhibited by Toupin and Gazis [7] through consideration of the nearest and next nearest interatomic interactions; they realized that the drawing in or pushing out the surface layer happens only in non-centrosymmetric materials. Later, Toupin, in a private communication with Mindlin [8], suggested that one can remove this restriction with the inclusion of the components of the second gradient of the strain tensor in the potential energy density function. Subsequently, Mindlin [8] proposed second strain gradient theory in which the strain energy density function depends on, not only the strain field and its first derivative, but also the second derivative of the strain field. Formulation within Mindlin’s second strain gradient theory gives rise to two surface parameters, namely, surface characteristic length and modulus of cohesion, enabling this theory to capture the surface effect of nano-structures on their mechanical properties [2,3,9]. Factually, consideration of the surface effect in nano-scale structures has been one of the most important stimuli in the development of higher order continuum theories. Recently, Shodja et al. [9] have proposed an atomistic model for the calculation of the additional constants for fcc materials in second strain gradient elasticity. Subsequently, they studied the surface effects on the behavior of nanosize Bernoulli-Euler beams. Zhang and Zhao [10] have noted that by detecting a very small shift in the resonance frequencies of a beam, it is possible to differentiate between the effects described by the nonlocality and those of the strain gradient and surface elasticities. By the employment of this dynamic method, they evaluate the unknown nonlocal parameters. Moreover, Ojaghnezhad and Shodja [3] employed a combined first principles and analytical approach for determination of the modulus of cohesion, surface energy, and the additional constants in second strain gradient elasticity. Later, Ojaghnezhad and Shodja [2] reformulated Gurtin and Murdoch [11] surface elasticity theory in the context of second strain gradient theory. In surface elasticity theory which has been formulated to capture the surface effects, the bulk material and surface layer are treated as two separate entities. Formulation in this framework entails the introduction of the notion of two surface parameters as surface residual stresses and surface elastic moduli tensor. The work of Ojaghnezhad and Shodja [2] leads to a linkage between the surface elastic parameters such as surface stress and surface elastic constants of surface elasticity theory and the elastic parameters stemming from second strain gradient elasticity. However, utilization of Mindlin’s second strain gradient theory to treat a vast variety of nano-scale problems involving various geometries in a mathematically rigorous manner requires the employment of suitable coordinate systems. Eringen [12] provides a simple but effective mathematical tool to transform any formulation written in the Cartesian coordinates to any curvilinear coordinate system. As an effort, Ji et al. [13] derived the general formulations of the simplified first strain gradient theory proposed by Zhou et al. [14] in the framework of orthogonal curvilinear coordinates. However, to date, in spite of the vast application of Mindlin’s second strain gradient theory in prediction of the mechanical behavior of nano-structures, its general formulation in orthogonal curvilinear coordinates is absent in the literature. In this paper, the methodology proposed by Eringen [12] is utilized to provide the formulation of Mindlin’s second strain gradient theory in any arbitrary orthogonal curvilinear coordinates. As it was alluded to, such formulations would be of great value for the treatment of nano-structures of various shapes where the surface effect is important. In addition, in order to illustrate the usefulness of the presented formulation in applications, relaxation of a spherical domain and a spherical cavity is examined based on Mindlin’s second strain gradient theory. Recently, solid and hollow nanospheres made of dielectric and precious metals such as gold and silver have absorbed great attention of the researchers due to their effective application in nanotechnology. Dielectric nanospheres are promising structures for light trapping in plannar thin-film solar cells [15]. Moreover, metal nanospheres due to their optical properties, have various technological applications such as surface plasmon resonance detection and imaging, surface-enhanced Raman scattering, and biomedical imaging and therapy. In the case of hollow gold nanospheres, the unique combination of small size, spherical shape, and strong tunable surface plasmon resonance is ideal for biomedical applications [16]. Hollow Pd spheres have been fabricated for usage as heterogeneous catalyst for suzuki coupling reactions [17]. Tunability of surface plasmon resonance by interior cavity size in Au hollow nanospheres has been examined by Liang et al. [18]. Au hollow nanosphere has also been used for drug delivery [19]. Metal nanoshell has found application in tumor therapy [20]. The present work focuses on the phenomenon of relaxation as an application of the current theoretical developments. In particular, we examine the relaxation of spherical domain as well as spherical cavity made of Ag, Au, and Pt based on Mindlin’s second strain gradient theory. The results are compared with the corresponding ones obtained from Gurtin–Murdoch surface elasticity as well as molecular dynamics simulation. 2. Second strain gradient theory in Cartesian coordinates Based on strain gradient theory formulated by Mindlin [8], strain energy density of a homogeneous and centrosymmetric material depends not only on the traditional infinitesimal strain,  ij , but also on its first and second spatial gradients,  ijk and  ijkl , respectively, as below

W =

1 1 1 C i j kl + Fi jklmn i j klmn + Gi jklmn i jk lmn + Ii jklmnpq i jkl mnpq + B◦i jkl i jkl . 2 i jkl 2 2

(2.1)

where the summation convention for repeated indices is employed and

1 2

i j = (ui, j + u j,i ),

(2.2a)

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

671

i jk = uk,i j ,

(2.2b)

i jkl = ul,i jk .

(2.2c)

In the above relations ui is the displacement component and “,” in subscript denotes the usual partial differentiation with respect to the Cartesian coordinates xi , i = 1, 2, 3. Based on the considered strain energy density, second-, third-, and ∂W ∂W ∂W forth-order stress tensors of any hyperelastic material are defined as τi j = , τi jk = , and τi jkl = . Thus,

∂i j

∂i jk

∂i jkl

τi j = Ci jkl kl + Fi jklmn klmn ,

(2.3a)

τi jk = Gi jklmn lmn ,

(2.3b)

τi jkl = Fpqi jkl  pq + Ii jklmntu mntu + B◦i jkl .

(2.3c)

b0 δ , where b0 is Mindlin’s modulus of cohesion and δi jkl = δi j δkl + δik δ jl + δ jk δil in which 3 i jkl δ ij is the Kronecker delta. Utilizing the above relations together with the symmetry considerations of the strain and stress tensors, it is inferred that the fourth, sixth, and eighth order tensors, respectively, Cijkl , Fijklmn , Gijklmn , and Iijklmntu have the following symmetry properties In the last relation, B◦i jkl =

Ci jkl = Ckli j = C jikl = Ci jlk ,

(2.4a)

Fi jklmn = Fjiklmn = Fi jlkmn = Fi jmlkn = Fi jkmln ,

(2.4b)

Gi jklmn = Glmni jk = G jiklmn = Gi jkmln ,

(2.4c)

Ii jklmnrs = Imnrsi jkl = I jiklmnrs = Ik jilmnrs = Ii jlkmnrs = Ii jklrnms = Ii jklnmrs = Ii jklmrns .

(2.4d)

Based on the above-mentioned symmetries and Eq. (2.2), the strain energy density function can be represented in the following form

W=

1 1 C ui, j uk,l + Fi jklmn ui, j un,klm + Gi jklmn uk,i j un,lm 2 i jkl 2 1 ◦ + Ii jklmnpq ul,i jk uq,mnp + Bi jkl ul,i jk . 2

(2.5)

For isotropic materials, the components of the fourth order tensor, Cijkl are written in terms of the usual Lamé constants

λ and μ as C1111 = λ + 2μ, C1122 = λ, and C1212 = μ. The not-mentioned nonzero components are obtained via cyclic permutation of indices. The nonzero components of the higher order elastic tensors, Fijklmn , Gijklmn , and Iijklmnpq for isotropic materials are related to Mindlin’s additional constants ai ’s, i = 1, . . . , 5, bi ’s, i = 1, . . . , 7, and ci ’s, i = 1, 2, 3 as displayed in Tables 1 and 2. The other nonzero components which are not displayed in the table can be obtained through the cyclic permutation of indices of the presented components. A material is referred to as “grade N” if the order of the highest position gradient in its energy density function expression is equal to N. Mindlin [8] showed that for a grade 3 material of volume V with boundary S, the stress-equation of motion in rectangular coordinate system has the following form

τip,i − τi jp,i j + τi jkp,i jk + f p = ρ u¨p ,

(2.6)

in which fp is the body force per unit volume and ρ is the mass density of the material. By substituting from Eqs. (2.3) and (2.2) into the stress-equation of motion, the displacement-equation of motion is derived as below

  ρ u¨ i = C jikl uk,l j + Fpq jkli + Fliq jkp − Gkli jqp u p,q jkl + I jklimnrs us,mnr jkl + fi .

(2.7)

For isotropic materials, the equation of motion is written in terms of Lamé constants and Mindlin’s additional parameters as below

      (λ + 2μ ) 1 − 211 ∇ 2 1 − 212 ∇ 2 u j, ji − μ 1 − 221 ∇ 2 1 − 222 ∇ 2 ei jk ekml ul,m j + fi = ρ u¨ i ,

(2.8)

672

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698 Table 1 The relations between the higher order tensors, Fijklmn and Gijklmn and Mindlin’s additional constants, ai ’s and ci ’s. c1 + c2 3 = c1

c1 + c3 3 c1 = 3 c3 = 6

F111122 = F111133 =

F111221 = F111331 =

F112222 = F113333

F112233 = F113322

c3 F121112 = F131113 = 2 c2 c3 F121121 = F131131 = + 3 6 c2 F121233 = 6 a2 G112233 = G113322 = 2 G221111 = G331111 = a2 + 2a3 a2 + 2a5 G112211 = G113311 = 2 G112332 = 2a3

F121332 = F122331

F111111 = c1 + c2 + c3 2a1 + a2 2 G112112 = G113113 = 2(a3 + a4 ) a1 + 2a4 + a5 G122122 = G133133 = 2 G123123 = a4 a1 G122133 = 2 G111122 = G111133 =

G111111 = a¯ a5 G123132 = 2

Table 2 The relations between the higher order tensor, Iijklmnpq and Mindlin’s additional constants, bi ’s. I11111111 = b¯ I11112222 = 2b1 I11122221 = 2b4 I11121112 = I11131113 = 2(b5 + b6 ) b3 + 2b5 3 b3 = I11131232 = 6 2b4 = I11133221 = 3 2 ( 2b2 + b3 + b4 ) = 9 b3 + 2b4 = I11311223 = 9 2 ( b1 + b2 ) = 9 2b1 + b3 = I11331221 = 9 2 ( b5 + 3b6 ) = I11321132 = 9 b3 + 4b7 = I11321231 = 18

I11121222 = I11131333 = I11121233 I11122331 I11211222 I11211332 I11221133 I11221331 I11231123 I11231231

2b1 + 2b2 + b3 3 2b1 + b3 + 2b4 + 2b5 I11111221 = I11111331 = 3 2b1 I11112233 = I11113322 = 3 b3 + 2b4 + 2b7 I11121121 = I11131131 = 3 b3 + 2b5 I11212331 = I11232333 = 9 2b5 I11121332 = I11131223 = 3 2 ( 2b2 + b3 + b5 + 3b6 + 2b7 ) I11211121 = I11311131 = 9 4b2 + b3 I11211233 = I11311322 = 18 2 ( b1 + b2 + b4 + b5 + 3b6 + b7 ) I11221122 = I11331133 = 9 2 ( b1 + b3 + 2b7 ) I11221221 = 9 2 ( b1 + b4 + b5 ) I11222332 = 9 2 ( b4 + b7 ) I11231132 = 9 b2 + 3b6 + b7 I12311231 = 9

I11111122 = I11111133 =

where eijk is the permutation tensor and

2(λ + 2μ ) 21 p = a¯ − 2c¯ ± [(a¯ − 2c¯ )2 − 4b¯ (λ + 2μ )] 2 ,

(2.9a)

2μ 22 p = a¯  − c3 ± [(a¯  − c3 )2 − 4b¯  μ] 2 ,

(2.9b)

1

1

for p = 1 and 2 pertinent to the positive and negative signs, respectively, and

a¯ = 2(a1 + a2 + a3 + a4 + a5 ),

(2.10a)

b¯ = 2(b1 + b2 + b3 + b4 + b5 + b6 + b7 ),

(2.10b)

c¯ = c1 + c2 + c3 ,

(2.10c)

a¯  = 2(a3 + a4 ),

(2.10d)

b¯  = 2(b5 + b6 ).

(2.10e)

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

673

11 , 12 , 21 , and 22 are the so-called “bulk characteristic lengths” which are related to Lamé constants and the additional parameters as given by Eq. (2.9). Moreover, Mindlin’s second strain gradient theory gives rise to another physically important length scale defined as

210 =

c¯

λ + 2μ

.

(2.11)

This length scale appears in the surface energy [2,8] and surface residual stress formula and thus, it is referred to as surface characteristic length [2]. Suppose that the outward unit normal at any point along S is defined by n(x). Then the generalized surface tractions, 1 2

3

t , t , and t on S are derived as 1

t i = n j (τ ji − τk ji,k + τkl ji,kl ) + Lk (n j τ jki − n j τm jki,m ) + Lk (L j (nm τm jki )) − Lt (nt nk nm n j n p,p τm jki − Lt (nk )nm n j τm jki ),

(2.12a)

t i = n j nk (τ jki − τl jki,l ) + nl Lk (n j τ jkli ) + Ll (nk n j τ jkli ),

(2.12b)

2

3

t i = n j nk nl τ jkli ,

(2.12c)

in which Li = ni n p,p − ∇i + ni n j ∇ j and ∇i = ∂ /∂ xi . 3. Second strain gradient theory in orthogonal curvilinear coordinates In this section, the stress-equation of equilibrium as well as the boundary conditions of second strain gradient theory described in the previous section is derived in the framework of orthogonal curvilinear coordinates. To this end, consider a set of orthogonal curvilinear coordinates xi , i = 1, 2, 3 with base vectors gi and metric tensor gi j = gi .g j [12]. The base vectors of the curvilinear coordinate system are obtained via gi = ∂ r/∂ xi where r = xk ik is the position vector in Cartesian coordinate system; ik , k = 1, 2, 3 are the Cartesian coordinate base vectors. So the square of the element of the arc length ds in terms of the curvilinear coordinates is written as below

ds2 = dr.dr = gkm dxk dxm .

(3.1)

In formulations within curvilinear coordinates, Einstein summation convention is applied for repeated indices on diagonal positions. Representation of the position vector r = xk ik in the above discussion is based on this convention. The unit base vectors of the curvilinear coordinates are obtained as

ei =

gi

|gi |

(no sum).

(3.2)

Now, consider the displacement vector u and express it via the curvilinear coordinates as

u = ui gi = u(i ) ei ,

(3.3)

where the physical components of u denoted by u(k) is obtained as

u(i ) = |gi |ui ,

(no sum).

(3.4)

According to Eringen [12], the passage from rectangular coordinates to curvilinear coordinates is made by replacing the usual partial differentiation by the covariant partial derivative which is indicated by the symbol “;”. Hence, the stress-equation of equilibrium in the curvilinear coordinates is displayed as follows

σ i p;i + f p = 0,

(3.5)

in which

σ i p = τ i p − τ i j p; j + τ i jk p; jk .

(3.6)

In the above relation, τ i p , τ ij p , and τ ijk p are the mixed components of the stress tensors of second-, third-, and fourthorder, respectively. According to Eringen [12], Ai... j k...l is called a mixed tensor if it changes under coordinate transformation through the following rule

Ai... j k...l (x ) = Am...n p...q (x )

∂ x i ∂ x j ∂ x p ∂ xq ... ... , ∂ xm ∂ xn ∂ x k ∂ x l

(3.7)

674

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

where x i denotes the component of the new coordinate system. Moreover, covariant partial differentiation of the second-, third-, and fourth-order mixed tensors is defined as below

 

σ i p; j = σ i p, j + σ t p

i jt

 

− σ it

t , jp

  i kt

σ i j p;k = σ i j p,k + σ t j p

 

σ

p;m

=σ

i jk

p,m

+σ

t jk

p

j kt

+ σ it p

 i jk

(3.8a)



 

t , kp

− σ i jt



i j + σ itk p km mt



(3.8b)

 +σ

i jt

p

k mt



 −σ

i jk

t



t , mp

where { jki } is the Christoffel symbol of the second kind defined as follows

  i jk

=

∂ 2 xn ∂ xi . ∂ x j ∂ xk ∂ xn

(3.9)

It can easily be shown that the physical components of the second-, third-, and fourth-order tensors are obtained as follows

τ (i ) ( j ) =

|gi | i τ , |g j | j

τ (i)( j ) (k) =

(no sum)

|gi ||g j | i j τ k, |gk |

(no sum)

|gi ||g j ||gk | i jk τ l. |gl |

τ (i)( j )(k) (l ) =

(3.10a)

(3.10b)

(no sum)

(3.10c)

According to the definition of covariant partial derivative (3.8), the stress equation of equilibrium (3.5) becomes

 

σ where

i

p,i

+σ

t

p

 

i t − σ it + f p = 0, ti ip

(3.11)

      i j t σ i p = τ i p − τ i j p, j − τ t j p − τ it p + τ i jt tj tj pj          i j k t + τ i jk p,k + τ t jk p + τ itk p + τ i jt p − τ i jk t kt

kt

tk

kp

,j

          i t j k n t jk n jk tnk t jn t jk + +τ p +τ p −τ n τ p,k + τ p tj kn nk nk kp           j i t k n itk ntk ink itn itk + +τ p +τ p −τ n τ p,k + τ p tj kn nk nk kp           t i j k n i jk n jk ink i jn i jk − +τ t +τ t −τ n . τ t,k + τ t jp

kn

nk

nk

kt

(3.12)

Subsequently, the equilibrium equation can be rewritten in terms of the physical components as below



|g p | ( i ) σ ( p) |gi |

where

σ

i

p



    |g p | i (i ) |gt | t + σ ( p) − σ (t ) + f ( p ) |g p | = 0, |gt | ti |gi | ip ,i (t )

   |g p | (i)( j ) |g p | ( i ) |g p | ( i ) = σ ( p) = τ ( p) − τ ( p) |gi | |gi | |gi ||g j | ,j j         |g p | (t )( j ) |g p | (i)(t ) i j |gt | (i)( j ) t − τ + τ − τ ( p) ( p) (t ) | gt ||g j | tj | gi ||gt | tj | gi ||g j | pj j,t    |g p | (i )( j )(t ) − τ ( p) |gi ||g j ||gt | , jt

(3.13)

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

     |g p | |g p | i j τ (t )( j )(k) ( p) + τ (i)(t )(k) ( p) | gt ||g j ||gk | kt | gi ||gt ||gk | kt ,j ,j j,t,k       |g p | k g | t | t + τ (i)( j )(t ) ( p) − τ (i)( j )(k) (t ) |gi ||g j ||gt | tk |gi ||g j ||gk | kp ,j ,j         |g p | |g p | i j + τ (t )( j )(k) ( p) + τ (i)(t )(k) ( p) |gt ||g j ||gk | tj | gi ||gt ||gk | tj ,k ,k     t |gt | − τ (i)( j )(k) (t ) |gi ||g j ||gk | jp ,k      |g p | t i + τ (n)( j )(k) | gn ||g j ||gk | kn t j ( p) n, j,k,t       |g p | |g p | j i k i + τ (t )(n)(k) ( p) + τ (t )( j )(n) ( p) |gt ||gn ||gk | nk t j |gt ||g j ||gn | nk t j       |g p | |gn | n i i j − τ (t )( j )(k) (n) + τ (n)(t )(k) ( p) |gt ||g j ||gk | kp t j |gn ||gt ||gk | kn t j       |g p | |g p | t j k j + τ (i)(n)(k) ( p) + τ (i)(t )(n) ( p) |gi ||gn ||gk | nk t j |gi ||gt ||gn | nk t j       |gn | n j |gt | i t − τ (i)(t )(k) (n) − τ (n)( j )(k) (t ) |gi ||gt ||gk | kp t j |gn ||g j ||gk | kn jp       |gt | j t |gt | k t − τ (i)(n)(k) (t ) − τ (i)( j )(n) (t ) |gi ||gn ||gk | nk jp |gi ||g j ||gn | nk jp    |gn | n t (i )( j )(k ) + τ . (n ) |gi ||g j ||gk | kt jp 

+

675



(3.14)

Moreover, the components of the second-, third-, and fourth-order strain tensors in curvilinear coordinates are written as

i j =

 





 

 1 i 1 i i m u ; j + g jm gin um ;n = u , j + ut + g jm gin um ,n + ut 2 2 tj nt

,

(3.15a)

          k k t k  k i j = uk ;i j = uk ,i + ut = uk ,i + ut + ut ,i + ur it

 −

uk ,t + ur

 l i jk = ul ;i jk  

 u ,t + u

r

 ul ,i + ur

+u

,j

t ir

 +

,j

,t

l

u ,t + u

r is

u

,t

+u

ur ,i + us

s

r ts

t − rj

;k

t ij

,k

    

ut ,r + us

t rs

r ij

l tk

    

 l

u ,r + u

 −

t ij

l tr



l − rj

l rt

l tr

  

ul ,t + ur

   r is

r



   r

 +

l − tj

  



   ur ,i + us

,j

 

l − tj

t ir

 +

l tr

t ir

  

ut ,i + ur

+

l ir

,i

r



  l

−

u

 

ut ,i + ur

+

−

,j

   t

+

l it



tj

(3.15b)



 

ul ,i + ut

ir

,j

t ij

l it

t

u ,i + u

=

k tr

  l

=

it

;j

  

ul ,r + us

s

l rs

r tj

t ik

    l rs

r it

t . kj

(3.15c)

676

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

Using Eqs. (3.4) and (3.10), the physical components of the second-, third-, and fourth-order strain tensors are derived.



(i )



1 ( j) = 2

          u(t ) |gi | i |g j | u ( j ) |g j | j |gi | u(i) + + + , |g j | |gi | , j |gi | |g j | ,i |gt | |g j | t j |gi | it t

(3.16a)

  (k)     (t )     u(t ) k |gk | u u k  (i)( j ) = + + |gi ||g j | |gk | ,i j | g | it | g | t j t t t ,j ,i  (k)   

 (r )       (r ) (k )

u |gk |

−

 (l ) (i)( j )(k) =

t ij

,t

+



u k k − tj |gr | tr

u t |gr | ir

t, r

t ij

(3.16b)

 (l )   (t )    (t )     |gl | u u l u l + + |gi ||g j ||gk | |gl | ,i jk | g | it | g | t j t t t , jk ,i ,k  (l )     (t )     (l )     (l )    u |gl |

− +





r, t

 +

 −

 −

u (r )

|gr | u (r ) |gr | u (r )

,t

t ij

+

,k

  

u (r ) t |gr | ir

l tj

    ,i

t rj

l − tk

l rj

t + ik

    ,t

    l

t

|gr | ,i rt k j     (s )

+

u |gt |

,i j



− ,k

  

u (l ) |gl | u (l )

u |gl |

,r

t ij

r tj

t − ik



t

u (r )

t kj

,it

   

u (r ) t |gr | ir

,j

    l

u l |gr | ir

+





u (s )

l tk

t

,t

t kj

    r

t

|gs | is r j r, s, t     (s )

−

u t |gs | rs

r ij

u l r − tk |gs | ts

l rj

t u l + ik |gs | rs

−

u (s ) r |gs | is

l rt

t u (s ) l + kj |gs | rs

r it

t kj

   

u |gl |

|gr | tr , j ik  (r )    

    r



,k

l − tk

|gl | ,r it k j     (s )

t − ik

+

r ij

    ,r

,t j

  

u (r ) l |gr | tr

   

u(t )

|gt |

l − tk

   

r tj

l tk

t ik

.

(3.16c)

Next, the traction boundary conditions in curvilinear coordinates are derived as below. To this end, define the reciprocal 3

base vectors gk such that gk .gl = δ k l where δ k l is the Kronecker delta. Moreover, gkm = gk .gm . Among the tractions, t has the simplest representations, and so we come up with its curvilinear representation first. 3

i

t = gis g jp gkq glr n p nq nr τ jkl s ,

(3.17)

and in terms of the physical components 3

(i )

 t n( p) n(q ) n(r ) ( j )(k )(l ) |gs | = gis g jp gkq glr τ . (s ) |gi | |g p | |gq | |gr | |g j ||gk ||gl |

(3.18)

j,k,l p,q,r,s

If orthogonal curvilinear coordinates are used, then it simplifies to 3

(i )

t

= n( j ) n(k ) n(l ) τ ( j )(k )(l ) (i ) .

(3.19)

The second traction boundary condition is similarly written as follows 2

m

t gmi = n p nq g p j gqk (τ k j i − τ l jk i;l ) + ns nr nx gsl grk gx j n p ;p τ jkl i − ns gsl (nr gr j τ jkl i );k + nq ns n p glq gsk (nx gx j τ jkl i );p + n p ;p nq nr ns gql grk gs j τ jkl i − (nr ns grk g js τ jkl i );l + n p nt g pl (nr ns grk gs j τ jkl i );t , and in terms of the physical components we have

(3.20)

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

2

677

(m )

t  n ( p ) n (q ) |gi | gmi = g p j gqk (τ (k )( j ) (i ) − τ (l )( j )(k ) (i );(l ) ) | g | | g | | g | | g m p q k ||g j | p, q m j, k, l

+

 n (s ) n (r ) n (x ) |gi | gsl grk gx j n( p) ;( p) τ ( j )(k )(l ) (i ) | g | | g | | g | | g | | s r x j gk ||gl | r, s, x j, k, l, p

 n (s ) |gi | − gsl (n(r ) gr j τ ( j )(k )(l ) (i ) );(k ) | g | | g | s j |gl ||gr | r, s j, k, l

+

 n (s ) n ( p ) n (q ) |gi ||g p | gsk gql (n(x ) gx j τ ( j )(k )(l ) (i ) );( p) | g | | g | | g | | g | |g j ||gk ||gl | s p q x q, r, s, x j, k, l, p

+



n ( p) ; ( p)

q, r, s j, k, l, p

−

n (s ) n (r ) n (q ) |gi | g g g τ ( j )(k )(l ) (i ) |gs | |gr | |gq | s j rk ql |g j ||gk ||gl |

 n (s ) n (r ) |gi | ( gs j grk τ ( j )(k )(l ) (i ) );(l ) | g | | g | | g s r j ||gk | r, s j, k, l

+

 n( p) n(t ) n (r ) n (s ) |gi ||gt | g pl ( grk gs j τ ( j )(k )(l ) (i ) );(t ) . | g | | g | | g | | g | | g p t r s j ||gk ||gl | p, q, r, s

(3.21)

j, k, l, t

In orthogonal curvilinear coordinates, it can further be simplified to 2

(i )

= n( j ) n(k ) (τ (k )( j ) (i ) − τ (l )( j )(k ) (i );(l ) ) + n(l ) n(k ) n( j ) n( p) ;( p) τ ( j )(k )(l ) (i )

t

− n(l ) (n( j ) τ ( j )(k )(l ) (i ) );(k ) + n(k ) n( p) n(l ) (n( j ) τ ( j )(k )(l ) (i ) );( p) + n( p) ;( p) n( j ) n(k ) n(l ) τ ( j )(k )(l ) (i ) − (n( j ) n(k ) τ ( j )(k )(l ) (i ) );(l ) + n(l ) n(t ) (n(k ) n( j ) τ ( j )(k )(l ) (i ) );(t ) ,

(3.22)

where for an arbitrary tensor of any order

A

(i )···( j )

( k ); ( m )

  |gi | · · · |g j | i··· j |gi | · · · |g j | i... j i t ... j = A k;m = A k,m + A k |gk ||gm | |gk ||gm | tm    

+ · · · + Ai...t k

|gi | · · · |g j | = |gk ||gm |

j t − Ai... j t tm km



|gk | A(i )···( j ) (k ) |gi | · · · |g j |  

 ,m

  |gk | i |gk | j + A(t )...( j ) (k ) + ··· + A(i )...(t ) (k ) |gt | · · · |g j | tm |gi | · · · |gt | tm   |gt | t − A(i )...( j ) (t ) . |gi | · · · |g j | km

(3.23)

Likewise, it can be shown that the first traction boundary condition in curvilinear coordinates which is more involved than the second and third traction types has the following representation 1

gmi t

m

= nt gt j (τ j i − τ k j i;k + τ kl j i;kl ) + n p ;p nt gtk (ns gs j τ jk i − ns gs j τ m jk i;m ) + nq gqk n p (nt gt j τ jk i − ns gs j τ m jk i;m );p − gm j nm ;k (τ jk i − τ p jk i;p ) − gt j nt (τ jk i;k − τ m jp i;mp ) + nt gtk ns ;s nl gl j n p ;p nn gnm τ m jk i − n p g pk ns ;s (nn gnm τ m jk i ); j + nr grk ns ;s nl gl j n p (nt gtm τ m jk i );p − (nt gt j n p ;p ns gsm τ m jk i );k + (nt gtm τ m jk i ); jk − (nt gt j n p (nr grm τ m jk i );p );k + nr grk ns (nl gl j n p ;p nx gxm τ m jk i );s − n p g pk ns (nq gqm τ m jk i ); js + nq gqk ns (nt gt j n p (nr grm τ m jk i );p );s − nt ns ;s n p ;t g pk nq gqm nr gr j τ m jk i

678

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

+ (n p ;h g pk ght nq nr gqm gr j τ m jk i );t − nt nr (n p ;t g pk nq nx gqm gx j τ m jk i );r + ns ;s n p nx gxm nt gt j nq ;p gqk τ m jk i − (nt n p nq gqm nx gx j ns ;p gsk τ m jk i );t + nx gtx ns (nt n p nq gmq nr gr j nz ;p gzk τ m jk i );s ,

(3.24)

which in terms of the physical components within the orthogonal curvilinear coordinate system has the following form 1

t

(i )

= n( j ) (τ ( j ) (i ) − τ (k )( j ) (i );(k ) + τ (k )(l )( j ) (i );(k )(l ) ) + n( p) ;( p) n(k ) n( j ) (τ ( j )(k ) (i ) − τ (m )( j )(k ) (i );(m ) ) + n(k ) n( p) (n( j ) τ ( j )(k ) (i ) − n( j ) τ (m )( j )(k ) (i );(m ) );( p) − n( j ) ;(k ) (τ ( j )(k ) (i ) − τ ( p)( j )(k ) (i );( p) ) − n( j ) (τ ( j )(k ) (i );(k ) − τ (m )( j )( p) (i );(m )( p) ) + n(k ) n(s ) ;(s ) n( j ) n( p) ;( p) n(m ) τ (m )( j )(k ) (i ) − n(k ) n(s ) ;(s ) (n(m ) τ (m )( j )(k ) (i ) );( j ) + n(k ) n(s ) ;(s ) n( j ) n( p) (n(m ) τ (m )( j )(k ) (i ) );( p) − (n( j ) n( p) ;( p) n(m ) τ (m )( j )(k ) (i ) );(k ) + (n(m ) τ (m )( j )(k ) (i ) );( j )(k ) − (n( j ) n( p) (n(m ) τ (m )( j )(k ) (i ) );( p) );(k ) + n(k ) n(s ) (n( j ) n( p) ;( p) n(m ) τ (m )( j )(k ) (i ) );(s ) − n(k ) n(s ) (n(m ) τ (m )( j )(k ) (i ) );( j )(s ) + n(k ) n(s ) (n( j ) n( p) (n(m ) τ (m )( j )(k ) (i ) );( p) );(s ) − n(t ) n(s ) ;(s ) n(k ) ;(t ) n(m ) n( j ) τ (m )( j )(k ) (i ) + (n(k ) ;(t ) n(m ) n( j ) τ (m )( j )(k ) (i ) );(t ) − n(t ) n(r ) (n(k ) ;(t ) n(m ) n( j ) τ (m )( j )(k ) (i ) );(r ) + n(s ) ;(s ) n( p) n(m ) n( j ) n(k ) ;( p) τ (m )( j )(k ) (i ) − (n(t ) n( p) n(m ) n( j ) n(k ) ;( p) τ (m )( j )(k ) (i ) );(t ) + n(t ) n(s ) (n(t ) n( p) n(m ) n( j ) n(k ) ;( p) τ (m )( j )(k ) (i ) );(s ) .

(3.25)

Finally, by substituting the components of the elastic tensors in Eq. (2.3) in terms of Lamé constants and Mindlin’s additional parameters, the constitutive relations for isotropic materials in the curvilinear coordinates are obtained as follows

τ ( p) (q) = λ (i) (i) δ ( p) (q) + 2μ ( p) (q) + c1  ( j ) (i)(i)( j ) δ ( p) (q) + c2  (i) ( p)(q)(i) +

c3 ( q ) ( (i)(i)( p) +  ( p) (i)(i)(q) ), 2

(3.26a)

a2 ( p ) ( (i)(i) δ (q) (r ) 2 + 2 (i ) (r )(i ) δ (q ) ( p) +  (q ) (i )(i ) δ ( p) (r ) ) + 2a3  (r ) (i )(i ) δ ( p) (q )

τ ( p)(q) (r ) = a1 ( (i) ( p)(i) δ (q) (r ) +  (i) (q)(i) δ ( p) (r ) ) +

+ 2a4  (r ) (q )( p) + a5 ( ( p) (r )(q ) +  (q ) (r )( p) ),

(3.26b)

2 3

τ ( p)(q)(r ) (s) = b1  ( j ) (i)(i)( j ) (δ ( p) (q) δ (r ) (s) + δ ( p) (r ) δ (q) (s) + δ (q) (r ) δ ( p) (s) ) 2 b2  (i ) ( j )(k )(i ) (δ ( j ) ( p) δ (k ) (q ) δ (r ) (s ) + δ ( j ) ( p) δ (k ) (r ) δ (q ) (s ) 3 1 + δ ( j ) (q ) δ (k ) (r ) δ ( p) (s ) ) + b3 ( (k ) (i )(i )( j ) (δ ( j ) ( p) δ (k ) (q ) δ (r ) (s ) 6 + δ ( j ) ( p) δ (k ) (r ) δ (q ) (s ) + δ ( j ) (q ) δ (k ) (r ) δ ( p) (s ) ) +  ( j ) (i )(i )(k ) (δ ( j ) ( p) δ (k ) (q ) δ (r ) (s ) +

+ δ ( j ) ( p) δ (k ) (r ) δ (q ) (s ) + δ ( j ) (q ) δ (k ) (r ) δ ( p) (s ) ) + 2 (i ) ( j )(s )(i ) (δ ( j ) ( p) δ (q ) (r ) 2 + δ ( j ) (q ) δ ( p) (r ) + δ ( j ) (r ) δ ( p) (q ) )) + b4  ( j ) (i )(i )(s ) (δ ( j ) ( p) δ (q ) (r ) 3 2 ( j) ( p) ( j) ( p) + δ (q ) δ (r ) + δ (r ) δ (q ) ) + b5  (s ) (i )(i )( j ) (δ ( j ) ( p) δ (q ) (r ) + δ ( j ) (q ) δ ( p) (r ) 3 2 ( j) ( p) (s ) + δ (r ) δ (q ) ) + 2b6  ( p)(q )(r ) + b7 ( ( p) (q )(r )(s ) +  (q ) (r )(s )( p) 3 1 (r ) (i ) ( p) (r ) +  (s )( p)(q ) ) + c1  (i ) (δ (q ) δ (s ) + δ ( p) (r ) δ (q ) (s ) + δ ( p) (s ) δ (q ) (r ) ) 3 1 (i ) (i ) + c2  ( j ) ( δ ( p ) δ ( j ) ( q ) δ ( r ) ( s ) + δ ( i ) ( p ) δ ( j ) ( r ) δ ( q ) ( s ) + δ ( i ) ( q ) δ ( j ) ( r ) δ ( p ) ( s ) ) 3 1 + c3  ( i ) ( s ) ( δ ( i ) ( p ) δ ( q ) ( r ) + δ ( i ) ( q ) δ ( p ) ( r ) + δ ( i ) ( r ) δ ( p ) ( q ) ) 3 1 + b0 ( δ ( p ) ( q ) δ ( r ) ( s ) + δ ( p ) ( r ) δ ( q ) ( s ) + δ ( p ) ( s ) δ ( q ) ( r ) ). 3

(3.26c)

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

679

Fig. 1. Geometrical representation of the Cartesian and spherical coordinate systems and the pertinent unit base vectors.

4. Second strain gradient theory in spherical coordinates In the spherical coordinate system shown in Fig. 1, the independent curvilinear variables x1 = r, x2 = θ , and x3 = φ are related to the Cartesian coordinates as

x1 = r sin φ cos θ , x2 = r sin φ sin θ , x3 = r cos φ . The corresponding base vectors are obtained via gi = ∂ r/∂ xi , i = 1, 2, 3. Thus, letting gr ≡ g1 , gθ ≡ g2 , and gφ ≡ g3 we obtain

gr = sin φ cos θ i + sin φ sin θ i + cos φ i , 1

2

3

gθ = −r sin φ sin θ i + r sin φ cos θ i , 1

2

gφ = r cos φ cos θ i + r cos φ sin θ i − r sin φ i . 1

2

3

(4.1)

The unit base vectors are obtained by the normalization of the above base vectors as

er = sin φ cos θ i + sin φ sin θ i + cos φ i , 1

2

3

eθ = − sin θ i + cos θ i , 1

2

eφ = cos φ cos θ i + cos φ sin θ i − sin φ i . 1

2

3

(4.2)

Since the base vectors given in Eq. (4.1) are orthogonal, then the nonzero diagonal components of the pertinent metric tensor and their corresponding reciprocal components are as below 2

grr = 1,

gθ θ = r 2 sin

grr = 1,

gθ θ =

φ,

1 r2

2

sin

φ

,

gφφ = r 2 ,

gφφ =

1 . r2

(4.3a)

(4.3b)

Moreover, the nonzero components of the Christoffel symbol of the second kind in spherical coordinates are as below



    r φ φ, = − sin φ cos φ , = −r, θθ φφ θθ             1 θ θ φ φ θ θ = = = = , = = cot φ . rθ rφ r θr φr θφ φθ r



2

= −r sin

(4.4)

680

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

Remark 1. In Section 2, the physical components of a quantity within a curvilinear coordinate system were indicated by embracing the pertinent superscripts and subscripts by parentheses. In what follows, for convenience we drop out the parentheses and, subsequently, move the superscript to the subscript. For example, the physical strain component  (i ) ( j ) with respect to the spherical coordinates will be presented as  rr ,  rθ ,  rφ ,  θ θ ,  φφ , and  φθ . Using Eq. (3.16), the strain field and its first and second gradients in the spherical coordinates, after some manipulations, are derived in terms of the components of the displacement field (ur , uθ , uφ ) as below

  1 1 rr = ur,r , rθ = ruθ ,r + ur,θ cscφ − uθ , rφ = ur,φ + ruφ ,r − uφ , 2r 2r   1 1 θ θ = ur + uθ ,θ cscφ + uφ cot φ , θ φ = u + uφ ,θ cscφ − uθ cot φ , r 2r θ ,φ  1 φφ = ur + uφ ,φ ,

(4.5a)

rrr = ur,rr , rrθ = uθ ,rr , rrφ = uφ ,rr ,     1 1 rθ r = 2 uθ − ur,θ cscφ + ur,rθ cscφ − uθ ,r , r r  1  1 rθ θ = ur,r + uθ ,rθ cscφ + uφ ,r cot φ − 2 ur + uθ ,θ cscφ + uφ cot φ , r r  1  1 rθ φ = uφ ,rθ cscφ − uθ ,r cot φ + 2 uθ cot φ − uφ ,θ cscφ , r r  1  1 1 1 rφ r = ur,rφ − uφ ,r + 2 uφ − ur,φ , rφθ = uθ ,rφ − 2 uθ ,φ , r r r r  1  1 rφφ = uφ ,rφ + ur,r − 2 ur + uφ ,φ , r r  1 1 2 θ θ r = 2 ur,θ θ csc φ − 2uθ ,θ cscφ + ur,φ cot φ − 2uφ cot φ − ur + ur,r , r r  1 1 2 θ θ θ = 2 uθ ,θ θ csc φ + 2uφ ,θ cot φ cscφ + 2ur,θ cscφ + uθ ,φ cot φ − uθ csc2 φ + uθ ,r , r r  1 1 2 2 θ θ φ = 2 uφ ,θ θ csc φ − 2uθ ,θ cot cscφ + uφ ,φ cot φ − uφ cot φ + uφ ,r , r r  1 θ φ r = 2 ur,θ φ cscφ − uφ ,θ cscφ − ur,θ cscφ cot φ − uθ ,φ + uθ cot φ , r  1 θ φθ = 2 ur,φ + uθ ,θ φ cscφ − uθ ,θ cot φ cscφ + uφ ,φ cot φ − uφ csc2 φ , r  1 θ φφ = 2 uφ ,θ φ cscφ − uφ ,θ cscφ cot φ + ur,θ cscφ − uθ ,φ cot φ + uθ cot2 φ , r  1 1 1 1 φφ r = ur,r + 2 ur,φφ − 2uφ ,φ − ur , φφθ = uθ ,r + 2 uθ ,φφ , r r r r  1 1 φφφ = uφ ,r + 2 uφ ,φφ + 2ur,φ − uφ ,

(4.5b)

r

r

r

rrrr = ur,rrr , rrrθ = uθ ,rrr , rrrφ = uφ ,rrr ,  1  2  2 rrθ r = 3 ur,θ cscφ − uθ + 2 uθ ,r − ur,rθ cscφ + ur,rrθ cscφ − uθ ,rr , r r r  2  2 rrθ θ = 3 ur + uθ ,θ cscφ + uφ cot φ − 2 ur,r + uθ ,rθ cscφ + uφ ,r cot φ r r  1 + u + uθ ,rrθ cscφ + uφ ,rr cot φ r r,rrθ  2  1 2 rrθ φ = 3 uφ ,θ cscφ − uθ cot φ − 2 uφ ,rθ cscφ − uθ ,r cot φ − uθ ,rr cot φ r r r

− uφ ,rrθ cscφ ,  2  1  2 ur,φ − uφ + 2 uφ ,r − ur,rφ + u − uφ ,rr , r r,rrφ r3 r 2 2 1 rrφθ = 3 uθ ,φ − 2 uθ ,rφ + uθ ,rrφ , r r r rrφ r =

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

rrφφ = rθ θ r =

rθ θ θ =

rθ θ φ =

rθ φ r =

rθ φθ =

rθ φφ =

rφφ r = rφφθ = rφφφ = θ θ θ r =

θ θ θ θ =

681

 2  1  2 ur + uφ ,φ − 2 ur,r + uφ ,rφ + ur,rr + uφ ,rrφ , 3 r r r  2 ur + 2uφ cot φ − ur,φ cot φ + 2uθ ,θ cscφ − ur,θ θ csc2 φ r3  1 1 + 2 ur,rθ θ csc2 φ − 2ur,r − 2uθ ,rθ cscφ − 2uφ ,r cot φ + ur,rφ cot φ + ur,rr , r r  2 2 2 uθ csc φ − uθ ,φ cot φ − 2ur,θ cscφ − 2uφ ,θ cscφ cot φ − uθ ,θ θ csc φ 3 r  1 1 + 2 uθ ,rθ θ csc2 φ − 2uθ ,r + 2ur,rθ cscφ + 2uφ ,rθ cot φ cscφ + uθ ,rφ cot φ − uθ ,r cot2 φ + uθ ,rr , r r  2 2 2 uφ cot φ − uφ ,φ cot φ + 2uθ ,θ cot φ cscφ − uφ ,θ θ csc φ r3  1 1 + 2 uφ ,rθ θ csc2 φ − 2uθ ,rθ cot φ cscφ + uφ ,rφ cot φ − uφ ,r cot2 φ − uφ ,r + uφ ,rr , r r  2 −ur,θ φ cscφ + uφ ,θ cscφ + ur,θ cot φ cscφ + uθ ,φ − uθ cot φ r3  1 + 2 uθ ,r cot φ − uθ ,rφ − ur,rθ cot φ cscφ − uφ ,rθ cscφ + ur,rθ φ cscφ , r  2 uφ csc2 φ − ur,φ − uφ ,φ cot φ + uθ ,θ cot φ cscφ − uθ ,θ φ cscφ 3 r  1 + 2 −uφ ,r csc2 φ + ur,rφ + uφ ,rφ cot φ − uθ ,rθ cot φ cscφ + uθ ,rθ φ cscφ , r  2 −uθ cot2 φ + uθ ,φ cot φ − ur,θ cscφ + uφ ,θ cot φ cscφ − uφ ,θ φ cscφ r3  1 + 2 uθ ,r cot2 φ − uθ ,rφ cot φ + ur,rθ cscφ − uφ ,rθ cot φ cscφ + uφ ,rθ φ cscφ , r  1  1 2 ur + 2uφ ,φ − ur,φφ − 2 2ur,r + 2uφ ,rφ − ur,rφφ + ur,rr , r r3 r  1 2 1 − 3 uθ ,φφ + 2 uθ ,rφφ − uθ ,r + uθ ,rr , r r r  2  1 2 uφ − 2ur,φ − uφ ,φφ − 2 2uφ ,r − 2ur,rφ − uφ ,rφφ + uφ ,rr , r r3 r 1 2 2 3uθ csc φ − 3uθ ,φ cot φ − 5ur,θ cscφ − 2ur,θ cot φ cscφ − 6uφ ,θ cot φ cscφ r3  3  + 3ur,θ φ cot φ cscφ − 3uθ ,θ θ csc2 φ + ur,θ θ θ csc3 φ + 2 ur,rθ cscφ − uθ ,r , r 1 3 − 3ur + (−5 cos φ + cos 3φ )csc3 φ uφ + 3ur,φ cot φ + 3uφ ,φ cot2 φ 3 4 r − 5uθ ,θ csc3 φ + 3uθ ,θ φ cot φ cscφ + 3ur,θ θ csc2 φ + 3uφ ,θ θ cot φ csc2 φ + uθ ,θ θ θ csc3 φ +



 3 ur,r + uφ ,r cot φ + uθ ,rθ cscφ , r2

1 3uθ cot φ csc2 φ − 3uθ ,φ cot2 φ − 2uφ ,θ cscφ − 5uφ ,θ cot2 φ cscφ r3

3  + 3uφ ,θ φ cot φ cscφ − 3uθ ,θ θ cot φ csc2 φ + uφ ,θ θ θ csc3 φ + 2 uφ ,rθ cscφ − uθ ,r cot φ , r    1  2 2 θ θ φ r = 3 2 + 3 cot φ uφ − 2 + cot φ ur,φ − 3uφ ,φ cot φ + ur,φφ cot φ r

θ θ θ φ =

+ 4uθ ,θ cot φ cscφ − 2uθ ,θ φ cscφ − 2ur,θ θ cot φ csc2 φ − uφ ,θ θ csc2 φ + ur,θ θ φ csc2 φ +

θ θ φθ =

 1 ur,rφ − uφ ,r , 2 r

1 2 cot φ csc2 φ uθ − 2uθ ,φ csc2 φ + uθ ,φφ cot φ − 2ur,θ cot φ cscφ r3 − 2uφ ,θ cscφ cot2 φ − 2uφ ,θ csc3 φ + 2cscφ ur,θ φ + 2 cot φ cscφ uφ ,θ φ



− 2 cot φ csc2 φ uθ ,θ θ + csc2 φ uθ ,θ θ φ +

1 uθ ,rφ , r2



682

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

θ θ φφ =

1 r3



− ur + 2uφ cot3 φ + ur,φ cot φ − uφ ,φ − 2uφ ,φ cot2 φ + uφ ,φφ cot φ

+ 4uθ ,θ cot2 φ cscφ − 2uθ ,θ φ cot φ cscφ + ur,θ θ csc2 φ



− 2uφ ,θ θ cot φ csc2 φ + uφ ,θ θ φ csc2 φ +

θ φφ r =

1 r3



 1 ur,r + uφ ,rφ , r2

− 2 cot2 φ uθ + 2uθ ,φ cot φ − uθ ,φφ − 2ur,θ cscφ + ur,θ cscφ cot2 φ

+ ur,θ csc3 φ + 2uφ ,θ cot φ cscφ − 2ur,θ φ cot φ cscφ − 2uφ ,θ φ cscφ + ur,θ φφ cscφ +

θ φφθ =

1 r3



 1 ur,rθ cscφ − uθ ,r , 2 r

− ur + cot φ (2csc2 φ − 1 )uφ − 2uφ ,φ csc2 φ + ur,φφ + uφ ,φφ cot φ



+ 2uθ ,θ cot2 φ cscφ − 2 cot φ cscφ uθ ,θ φ + uθ ,θ φφ cscφ +



1 ur,r + uφ ,r cot φ r2

+ uθ ,rθ cscφ ,

θ φφφ =

1 r3



− 2 cot3 φ uθ + 2uθ ,φ cot2 φ − uθ ,φφ cot φ − 2ur,θ cot φ cscφ − 2uφ ,θ cscφ

+ uφ ,θ cscφ cot2 φ + uφ ,θ csc3 φ + 2ur,θ φ cscφ − 2uφ ,θ φ cot φ cscφ



 1 uφ ,rθ cscφ − uθ ,r cot φ , r2  3  1 φφφ r = 3 3uφ − 5ur,φ − 3uφ ,φφ + ur,φφφ + 2 ur,rφ − uφ ,r , r r  3 1 φφφθ = 3 uθ ,φφφ − 2uθ ,φ + 2 uθ ,rφ , r r  3  1 φφφφ = 3 −3ur − 5uφ ,φ + 3ur,φφ + uφ ,φφφ + 2 ur,r + uφ ,rφ . r r + uφ ,θ φφ cscφ +

(4.5c)

Finally, the set of equilibrium equations in spherical coordinate system is derived as follows

τ¯rr,r +

 1 1 1 τ¯θ r,θ + τ¯φ r,φ + 2τ¯rr + τ¯φ r cot φ − τ¯θ θ − τ¯φφ + fr = 0, r sin φ r r

(4.6a)

τ¯rθ ,r +

 1 1 1 τ¯ + τ¯φθ ,φ + 2τ¯rθ + τ¯θ r + (τ¯θ φ + τ¯φθ ) cot φ + fθ = 0, r sin φ θ θ ,θ r r

(4.6b)

τ¯rφ ,r +

 1 1 1 τ¯θ φ ,θ + τ¯φφ ,φ + 2τ¯rφ + τ¯φ r + (τ¯φφ − τ¯θ θ ) cot φ + fφ = 0, r sin φ r r

(4.6c)

where the expressions for τ¯i j , i, j = r, θ , φ are given in Appendix A. 5. Application to problems involving spherical symmetry As a three-dimensional useful application, consider a spherical surface of radius a bounding an isotropic body. In the absence of any type of external loading and under the assumption of central symmetry, it is expected that due to the effect of surface, the components of the displacement field along θ and φ be zero, uθ = uφ = 0 and hence

u = ( ur ( r ), 0, 0 ).

(5.1)

Based on the formulations given in Section 4 and using the relations (4.6), the only nontrivial equilibrium equation governing the body has the following form 2 2  r 2 211 212 u + 6r211 212 u − (r 2 (211 + 212 ) + 6211 212 )u r r r − 4r (11 + 12 )ur

+ (4(211 + 212 ) + r 2 )ur + 2rur − 2ur = 0.

(5.2)

By defining the new parameters p1 and p2 as

p1 =

211 + 212 211 212

,

p2 =

1 , 211 212

the equilibrium equation may be rewritten as 5  4 6  5  4 r 6 u + 6r 5 u − (6r 4 + p1 r 6 )u r r r − 4r p1 ur + ( 4 p1 r + p2 r )ur + 2 p2 r ur − 2 p2 r ur = 0.

(5.3)

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

683

In the above equation, r = 0 is a regular singular point. Its solution may be obtained by employing the Frobenius series ∞ 

ur =

an r n+s .

(5.5)

n=0

Its substitution into Eq. (5.4) yields ∞ 

(n + s )(n + s − 1 )(n + s − 2 )(n + s − 3 )(n + s − 5 )(n + s + 2 )an rn+s

n=0

− p1 + p2

∞  n=2 ∞ 

(n + s − 2 )(n + s − 3 )(n + s − 5 )(n + s )an−2 rn+s (n + s − 5 )(n + s − 2 )an−4 rn+s = 0.

(5.6)

n=4

The indicial polynomial, P(s) corresponding to n = 0 becomes

P (s ) = s(s − 1 )(s − 2 )(s − 3 )(s − 5 )(s + 2 ).

(5.7)

By the assumption a0 = 0, s is required to be a root of the indicial polynomial. It is observed that P(s) has six integer roots as s = −2, 0, 1, 2, 3, and 5. In Eq. (5.6), by equating the coefficients of r n+s for n ≥ 1, the following relations are resulted

n=1:

a1 (s + 1 )s(s − 1 )(s − 2 )(s − 4 )(s + 3 ) = 0,

(5.8a)

n=2:

s(s − 1 )(s − 3 )(s + 2 )[(s + 1 )(s + 4 )a2 − p1 a0 ] = 0,

(5.8b)

n=3:

s(s + 1 )(s − 2 )(s + 3 )[(s + 2 )(s + 5 )a3 − p1 a1 ] = 0,

(5.8c)

an P (n + s ) − p1 an−2 (n + s − 2 )(n + s − 3 )(n + s − 5 )(n + s )

n≥4:

+ p2 an−4 (n + s − 5 )(n + s − 2 ) = 0.

(5.8d)

A close scrutiny of the above relations reveals that by considering the roots, s = 2, −2, 0 the complete list of the indep1 a0 , a3 = 0, 3∗6 and pendent solutions will be obtained. First of all by considering s = 2 relations (5.8a)–(5.8d) result in a1 = 0, a2 =

a2n =

a2n+1

−2a0 p2 (n + 1 )(n + 2 )(2n + 1 )!

1680a3 n + 2 =− p 2 ( 2 n + 5 )!

n+1 2  k=1

n+2 2 



k=1

 (−1 )

k n+1−2k k p1 p2



n≥2

(5.9a)



n−k k−1



60a1 (n + 2 )  n−1−k + (−1 )k pn1−2k pk2 , ( 2 n + 5 )! k−1 n/2



n+1−k , k−1

(−1 )k pn1+2−2k pk2

n ≥ 2.

(5.9b)

k=1

Subsequently, employing the relation (5.5), the corresponding solution reduces to



ur = −30(211 + 212 )a1 840a3 + p2 2 a0 − p2





6 6 ur1 − 4 11 4 ur2 + 4 12 4 ur3 3 11 − 12 11 − 12

4 4 ur1 − 2 11 2 ur2 + 2 12 2 ur3 3 11 − 12 11 − 12







411 412 411 412 4ur4 + 2 ur − ur − ur + ur , 11 − 212 5 211 − 212 2 211 − 212 6 211 − 212 3

(5.10)

where a0 , a1 , and a3 are arbitrary and

ur1 = r,

(5.11a)

684

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

r 11

cosh ur2 =

r 12

cosh ur3 =

ur4 =

−

r

−

r

11 sinh

r 11

r2 12 sinh r2

r 12

,

(5.11b)

,

(5.11c)

1 , r2

(5.11d)

r  (11 + r )e 11 ur5 = , r2 −

(5.11e)

r  (12 + r )e 12 ur6 = . r2 −

(5.11f)

Likewise, for s = −2 it can be shown that



ur = a0 +

411 412 411 412 −3ur4 + 2 ur − ur − ur + ur 11 − 212 5 211 − 212 2 211 − 212 6 211 − 212 3



6 a3 2a2 11 12 (12 ur5 − 12 ur2 − 11 ur6 + 11 ur3 ) − (ur − ur3 ), 211 − 212 2 p2 (211 − 212 ) 2

(5.12)

for arbitrary a0 , a2 , and a3 . Finally, if one assumes s = 0, the series solution will collapse to the following form





ur = a0 2 (

211

+

 + a2 +

212

3  2 )ur4 − 2  ur − 311 ur2 − 312 ur6 + 312 ur3 11 − 212 11 5

411 412 411 412 4ur4 + 2 ur − ur − ur + ur 11 − 212 5 211 − 212 2 211 − 212 6 211 − 212 3



 2  3a1  4 60a3 11 ur2 − 412 ur3 + 11 ur2 − 212 ur3 , 2 2 2 − 12 2 p2 (11 − 12 )

(5.13)

211

with arbitrary a0 , a1 , a2 and a3 . In view of Eqs. (5.10), (5.12), and (5.13) it is concluded that the solution to Eq. (5.4) associated with spherical geometry problem having central symmetry can be represented as

u r = L{ u r 1 , u r 2 , u r 3 , u r 4 , u r 5 , u r 6 } .

(5.14)

Next, the traction boundary conditions associated with the proposed spherical problem are extracted from the general treatment given in the earlier sections in curvilinear coordinates. Recall that the spherical body is traction free and deforms just under the surface effect. The nontrivial components of the tractions in second strain gradient theory for this special problem which are derived using the relations (4.5), (3.26), (3.25), (3.22), and (3.19) must be equalized to zero on the spherical free surface (r = a ) as below 1 tr

2b0 =− 2 + r



+

 + +

2 tr





4(a4 + a5 + a¯ − 2c2 − 2c3 ) 8(b¯ − 2b1 + 4b6 + 4b7 ) 2λ − − ur r r3 r5



4(a4 + a5 + a¯ + c1 ) − 10c¯ 8(b¯ − 2b1 + 4b6 + 4b7 ) + ur,r λ + 2μ + r2 r4



4c¯ − 2a¯ + 2c1 8 ( b2 + b3 + b4 + b5 ) + ur,rr + r r3

4b¯ ur,rrrr + b¯ ur,rrrrr = 0, r

2b0 = + r



on





2c¯ − a¯ −

10b¯ ur,rrr r2

r = a,

(5.15a)



4(a1 + a2 + a3 ) + 4c1 + 2c¯ 8(b¯ − 2b1 + 4b6 + 4b7 ) + ur r2 r4

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

 +

 + −

4(a1 + a2 + a3 ) 8(b¯ − 2b1 + 4b6 + 4b7 ) − ur,r r r3



a¯ − c¯ +

12b¯ − 8(b2 + b3 + b4 + b5 ) ur,rr r2

2b¯ ur,rrr − b¯ ur,rrrr = 0, r



3

t r = b0 +

 +

685



on

r = a,



2 c1 8 ( b2 + b3 + b4 + b5 ) + ur + r r3

(5.15b)





c¯ −



2b¯ + 4b1 − 4b6 − 4b7 ur,rr + b¯ ur,rrr = 0 r

8 ( b2 + b3 + b4 + b5 ) ur,r r2 on

r = a.

(5.15c)

For illustration, the formulations for spherical geometries are specialized to an isotropic solid nanosphere in Section 5.1 and isotropic infinite body with spherical nanocavity in Section 5.2. The corresponding numerical examples will be given in Section 7. 5.1. Isotropic solid nanosphere For an elastic isotropic spherical solid considered under the surface effect, the solution must be bounded at its center, Thus, from the solutions presented by Eq. (5.11), the solution of the problem of interest may be written as

ur = A1 ur1 + A2 ur2 + A3 ur3 ⎛ r r cosh 11 sinh   11 11 = A1 r + A2 ⎝ − r r2

⎞

⎛

⎠ + A3 ⎝

cosh

r 12

r

−

12 sinh r2

r 12

⎞ ⎠,

(5.16)

The unknown coefficients A1 through A3 are determined using the boundary conditions given by Eq. (5.15). 5.2. Infinite isotropic domain with spherical nanocavity For an infinite isotropic domain containing a spherical void under the surface effect, ur1 , ur2 , and ur3 given in Eq. (5.11) are not suitable since they do not diminish as r → ∞, and thus the displacement field has the following form

ur = A1 ur4 + A2 ur5 + A3 ur6 r r − −   11 A1  + r e  + r e ( 11 ) ( 12 ) 12 = 2 + A2 + A3 . r r2 r2

(5.17)

Again, in a similar manner to the previous case, the unknown coefficients A1 through A3 are determined using the boundary conditions given by Eq. (5.15). 6. Comparison to Gurtin–Murdoch surface elasticity solution One can obtain nontrivial solutions for two problems of isotropic spherical solid and infinite isotropic domain containing a spherical void under the surface effect in the framework of Gurtin–Murdoch surface elasticity [11], as well. According to Gurtin and Murdoch [11], in the absence of body forces, the system of governing equilibrium equations in the bulk volume V of the solid can be simplified as below



τ pq,p = 0, τ pq = λii δ pq + 2μ pq ,

(6.1)

where  ij follows definition given by Eq. (2.2a). In Gurtin-Murdoch theory, these equations are coupled with the following governing equilibrium equations on the surface of the volume, ∂ V

⎧ ⎨I jM SiM, j = n p τ pi ,

 1 PMi ui, j I jN + PM j ui, j IiN , 2 = σ0 IiM + λ0 ENN IiM + 2μ0 IiN ENM + σ0 ui, j I jM .

EMN =

⎩

SiM

(6.2)

In the above equations, upper-case indices belong to the two-dimensional subspace of the three-dimensional Euclidean space. Consider a three-dimensional vector space V and its two-dimensional subspace Tx through which the structure of a

686

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

surface ∫ at each x ∈ ∫ is defined. According to the definition given by Gurtin and Murdoch [11], IiM is an inclusion map that linearly transforms any vector in the two-dimensional subspace (Tx ) to its corresponding vector in the three-dimensional space (V ), while PMi is the perpendicular projection from V onto Tx . In the above relations, EMN is the tangential surface strain tensor, SiM is the first Piola-Kirchhoff surface stress tensor, ni is the outward unit normal vector of the surface ∫, and σ 0 , λ0 , and μ0 are, respectively, residual surface tension and surface Lamé moduli for the isotropic surface ∫. It may be noteworthy to mention that I jM SiM, j in the first of Eq. (6.2) represents the surface divergence of the surface stress tensor SiM as defined by Gurtin and Murdoch [11,21]. For the problems involving spherical symmetry, as discussed in Section 5, the displacement field has the form (5.1). The displacement form of the governing equilibrium equation for the bulk is derived using Eq. (6.1)

r 2 ur,rr + 2rur,r − 2ur = 0,

(6.3)

with the general solution

ur = A r +

B . r2

(6.4)

The unknown constants A and B will be obtained using the surface boundary conditions given by Eq. (6.2). 6.1. Isotropic solid nanosphere In this case, since r = 0 is a field point it requires that B = 0, and the unknown A is determined via the first relation of Eq. (6.2) on r = a. Thus, the displacement field will take the following form

ur =

−2σ0 r . 2(2λ0 + 2μ0 + σ0 ) + (3λ + 2μ )a

(6.5)

It is interesting to note that comparison of the above solution with that of Mindlin’s second strain gradient theory given by (5.16) reveals that the two solutions share the linear term, but Mindlin’s theory gives rise to two additional terms. The effect of these additional terms will be numerically demonstrated in Section 7. 6.2. Infinite isotropic domain with spherical nanocavity Since the displacement field induced by the surface due to the spherical nanocavity should diminish at infinity, then A = 0. By employing the relations (6.2) on r = a, the displacement field will have the following form

ur =

−σ0 a3 . (2μa + 2λ0 + 2μ0 + σ0 )r2

(6.6)

As it is seen, in this case surface elasticity theory recovers only the first term of the solution obtained via Mindlin’s second strain gradient theory. In fact, the Mindlin’s solution (5.17) contains two additional terms which are absent in (6.6). 7. Numerical results For the illustration of the current theoretical developments, surface relaxation of spherical domains are considered. More specifically, the displacement field of a nano-spherical medium as well as the displacement field within an infinite domain containing a spherical nanocavity is examined. To the end of comparison, the numerical results are given in the framework of both Mindlin’s second strain gradient elasticity and Gurtin-Murdoch surface elasticity. The results are presented for some face-centered cubic crystals (fcc) as Ag, Au, and Pt. First, a brief explanation for the determination of the numerical values of the material parameters in second strain gradient elasticity and Gurtin–Murdoch surface elasticity is given in Section 7.1. 7.1. Evaluation of Mindlin’s material parameters via lattice dynamics and ab-initio calculations In order to present the numeric solution of problems involving spherical symmetry through second strain gradient elasticity, one should first determine the pertinent numerical values of the components of the fourth, sixth, and eighth order elastic moduli tensors of the crystals of interest. To this end, a short introduction to the atomistic description of materials via lattice dynamics is given here. Consider the bulk of a centrosymmetric crystal with perfect lattice of infinite extension in space and denote the position of an arbitrary primitive unit cell of volume v within by the vector x. Suppose that the distance between the α th primitive unit cell from the reference unit cell at x is indicated as Rα . Moreover, let Kiαj present the atomic force constant between the unit cells with location vectors x and x + Rα . For any perturbation of the atomistic configuration from the equilibrium, the potential energy density function pertinent to the one-atom unit cell at x to within a harmonic approximation may be expressed as

W =−

  1  α Ki j (ui (x + Rα ) − ui (x ) ) u j (x + Rα ) − u j (x ) . 4v α

(7.1)

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

687

Fig. 2. Variation of the displacement field ur in A˚ versus r/a in the (a) Ag, (b) Au, and (c) Pt nanospheres of radius a according to (1) Mindlin’s second a . a0

strain gradient theory and (2) Gurtin–Murdoch surface elasticity for different ratios of α =

By writing Taylor’s expansion of u(x + Rα ) about x and based on the fact that for centrosymmetric crystals the odd-ranked elastic moduli tensors are equal to zero, the potential energy density function may be written as

W=

1 1 C˜i jmn ui,m u j,n + C˜i jmnpq ui,m u j,npq + C˜˜i jmnpq ui,mn u j,pq 2 2 1 ˜ + Ci jmnt pqr ui,mnt u j,pqr . 2

(7.2)

The coefficients C˜ appearing in the above relation depend on the atomic force constants and the equilibrium positions of the atoms as follows

C˜i jmn = −

1  α K Rα Rα , 2v α i j m n

C˜i jmnpq = −

1  α K Rα Rα Rα Rα , 12v α i j m n p q

(7.3a)

(7.3b)

688

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

Fig. 3. Variation of the displacement field ur in A˚ versus r/a in the (a) Ag, (b) Au, and (c) Pt infinite domain with spherical cavity of radius a according to a . a0

(1) Mindlin’s second strain gradient theory and (2) Gurtin–Murdoch surface elasticity for different ratios of α =

C˜i jmnpqrs = −

1  α K Rα Rα Rα Rα Rα Rα , 72v α i j m n p q r s

(7.3c)

3 C˜ . Employing the Hamilton’s principle as described by Ojaghnezhad and Shodja [2] and Shodja et al. 2 i jmnpq [22], the equations of motion for centrosymmetric crystals are obtained as and C˜˜i jmnpq =

ρ u¨ i = C˜i jmn u j,mn + C˜i jmnpq u j,mnpq + C˜i jmnpqrs u j,mnpqrs ,

(7.4)

in which ρ is the ratio of the mass of the atom in one primitive unit cell to its volume, and “,” in the subscript denotes differentiation with respect to x. The atomic force constants, Kiαj are equivalent to the components of the Hessian matrix which are in turn equal to the value of the second derivative of the total potential energy with respect to the corresponding atomic positions at the

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

689

a

a ) with the normalized radius of the sphere ( ) obtained via Mindlin’s second strain a0 a0 gradient theory (MSGT), Gurtin–Murdoch surface elasticity (GMSE), and LAMMPS simulation for (a) Ag, (b) Au, and (c) Pt. Fig. 4. Variation of the normalized change in the sphere radius (

equilibrium. The Hessian matrix is obtained from the first principles density functional theory (DFT) and, subsequently, the fourth, sixth, and eighth order constants given by relations (7.3) are evaluated. From comparison of Eqs. (2.7) and (7.4) in the absence of body forces fi the following relations are obtained

Cim jn = C˜i jmn + C˜m jin − C˜mi jn ,

C˜i jklmn =

(7.5a)

1 F + Fjlmnki + Fjmknli + Fjnmkli + Flikmn j + Fkilmn j + Fmikln j 4 jkmnli  1 + Fnikml j − Gnlimk j + Gnkiml j + Gklimn j + Gnmilk j + Gmlink j  6 + Gmkinl j ,

C˜i jklmnpq =

(7.5b)

1  I + Iqmliknp j + Iqnlimkp j + Iqplimnk j + Iqkmilnp j + Iqkniml p j 20 qklimnp j + Iqkpimnl j + Imkliqnp j + Inklimqp j + I pklimnq j + Iqmnikl p j + Iqmpikln j + Iqnpimkl j + Imkniql p j + Imkpiqnl j + Inkpimql j + Imnliqkp j + Impliqnk j



+ Inplimqk j + Imnpiqkl j .

(7.5c)

Using the above relations, a set of equations for a¯ − 2c¯, a¯  − c3 , b¯ , b¯  , λ, and μ, pertinent to isotropic materials, in terms of the components of the tensors C˜ is obtained. Subsequently, the bulk characteristic lengths are readily available via Eq. (2.9). To obtain the other additional parameters of Mindlin’s theory, one may equalize the higher order terms of the strain energy density functions pertinent to the continuum model (2.5) and the pertinent lattice dynamics formulation (7.2) as below

C˜˜kni jlm uk,i j un,lm = Gi jklmn uk,i j un,lm ,

(7.6a)

690

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

a

a ) with the normalized radius of the sphere ( ) obtained via a0 a0 Mindlin’s second strain gradient theory (MSGT) and Gurtin–Murdoch surface elasticity (GMSE) for (a) Ag, (b) Au, and (c) Pt. Fig. 5. Variation of the normalized change in the embedded spherical void radius (

Table 3 ˚ and modulus of Lamé constants in units of eV/A˚ 3 , bulk and surface characteristic lengths in units of A, cohesion in units of eV/A˚ for Ag, Au, and Pt. element

λ (eV/A˚ 3 )

μ (eV/A˚ 3 )

˚ 11 , 12 (A)

˚ 21 , 22 (A)

˚ 10 (A)

˚ b0 (eV/A)

Ag Au Pt

0.56 1.08 1.55

0.24 0.26 0.60

0.91 ± 1.03i 0.56 ± 0.63i 0.81 ± 0.91i

1.37 ± 1, 66i 0.69 ± 0.69i 1.44 ± 1.45i

2.16i 0.259+0.420i 0.271+0.936i

-1.87 -0.157+0.314i -0.984+0.622i

C˜in jklm ui, j un,klm = Fi jklmn ui, j un,klm ,

(7.6b)

C˜lqi jkmnp ul,i jk uq,mnp = Ii jklmnpq ul,i jk uq,mnp .

(7.6c)

Using the above equalities, one may obtain the following relations between the additional parameters

a1 = a2 = a5 ,

2a3 = a4 ,

(7.7a)

2 c1 = c2 = c3 ,

4b1 = 2b2 = b3 = 4b4 = 2b7 ,

(7.7b)

b5 =

3 b6 . 2

(7.7c)

Based on the numerical values of Lamé constants, bulk and surface characteristic lengths and modulus of cohesion for Ag, Au, and Pt given by Ojaghnezhad and Shodja [2] and summarized in Table 3, one may determine all the material parameters in Mindlin’s second strain gradient theory. Moreover, Ojaghnezhad and Shodja [2] have also provided the surface

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

691

Fig. 6. Relaxation phenomenon observed for the Ag spherical domains of radii (a) a = 0.71a0 , (b) a = a0 , and (c) a = 2a0 via continuum theories of interest (MSGT and GMSE) as well as some EAMs.

residual stress and surface elastic constants for Ag, Au, and Pt within Gurtin-Murdoch surface elasticity. For convenience, the numerical values of these are displayed in Table 4.

7.2. Descriptive examples Based on the numerical values of the material parameters and characteristic lengths as determined in the previous section, one may evaluate the displacement and stress field through the domains of the spherical symmetry in the framework

692

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

Fig. 7. Variation of the stress field components τ rr and τθ θ = τφφ in eV/A˚ 3 versus r/a in the (a) Ag, (b) Au, and (c) Pt solid nanosphere of radius a according a . to Mindlin’s second strain gradient theory for different ratios of α = a0

Table 4 Surface residual stress and surface elastic constants in units of eV/A˚ 2 for Ag, Au, and Pt. element

σ0

λ0

μ0

Ag Au Pt

0.088 0.073 0.160

-0.047 -0.028 -0.083

-0.044 -0.036 -0.081

of Mindlin’s second strain gradient theory. Let a0 denote the lattice constant for the element under consideration, then the a normalized parameter, α = provides a sense on the size of the spherical domain as compared to the lattice constant a0 of the element. Exploiting the given numerical data for the material constants, the variation of the displacement field ur in units of A˚ versus the normalized variable r/a for nanosphere and spherical nanocavity of radius a pertinent to various

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

693

a Fig. 8. Stress components τrr = τθ θ = τφφ in eV/A˚ 3 versus via Gurtin–Murdoch surface elasticity for solid nanospheres made of (a) Ag, (b) Au, and (c) a0 Pt.

a is plotted in Figs. 2 and 3, respectively. These plots are given for fcc metals of (a) Ag, (b) Au, and (c) Pt a0 according to (1) Mindlin’s second strain gradient elasticity and (2) Gurtin–Murdoch surface elasticity. Since Mindlin’s bulk characteristic lengths are complex numbers, the displacement fields pertinent to Mindlin’s second strain gradient theory r represent oscillatory behavior with increasing (Figs. 2 and 3). Moreover, it is observed that the Gurtin–Murdoch solutions a for nanosphere and spherical nanocavity are less sensitive to the sphere radius in comparison to Mindlin’s solutions. Additionally, Mindlin’s second strain gradient theory provides values of ur on the boundary of the nanosphere much larger than Gurtin–Murdoch surface elasticity.   a

a The normalized change of the radius versus the normalized radius of the sphere is also plotted for the fcc a0 a0 crystals of (a) Ag, (b) Au, and (c) Pt in Fig. 4(a)–(c) according to both Mindlin’s second strain gradient elasticity (MSGE) and Gurtin–Murdoch surface elasticity (GMSE). For the sake of comparison, the relaxation of the nanosphere is also calculated by simulating the spherical domain via the molecular dynamics package LAMMPS at absolute temperature using EmbeddedAtom-Method (EAM) functions reported by Foiles et al. [23]. It is observed that the result obtained from Gurtin–Murdoch surface elasticity is approximately size independent while theresults  of Mindlin’s strain gradient theory and LAMMPS are

a size-dependent. Likewise, the normalized change of the radius of an embedded nano-sized spherical cavity is plotted a0 a

versus the normalized cavity radius for Ag, Au, and Pt in Fig. 5(a)–(c), respectively, using Mindlin’s strain gradient a0 theory and Gurtin–Murdoch surface elasticity. As it is observed the phenomenon of relaxation captured within MSGT is remarkably affected by the size of the spherical cavity, whereas GMSE remains nearly size independent. In order to compare the results on relaxation of a spherical solid predicted by the current continuum theories (MSGT and GMSE) with those of atomistic simulations, the atomic displacements calculated for Ag spherical domains via LAMMPS using different EAMs [23–26] as well as the corresponding results of the thories of interest are inserted in a common picture, Fig. 6. In Fig. 6(a)–(c), three different sizes of solid nanospheres with ratios α = 0.71, 1, and 2 are considered, respectively. In Fig. 6(a)–(b), it was feasible to indicate the atoms in a common radial distance from the center atom by letters. The computed radial displacement for each atom is also given in A˚ in this figure. Under the surface effect in the above-discussed problems involving spherical symmetry, the non-trivial displacement field implies non-trivial stress field with non-zero components τ rr (r) and τθ θ (r ) = τφφ (r ) in the unrelaxed configuration. values of α =

694

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

Fig. 9. Variation of the stress field components τ rr and τθ θ = τφφ in eV/A˚ 3 versus r/a in the (a) Ag, (b) Au, and (c) Pt domain with spherical cavity of a . radius a according to Mindlin’s second strain gradient theory for different ratios of α = a0

r Fig. 7(a)–(c) represent the variation of τ rr and τθ θ = τφφ in eV/A˚ 3 versus for some different values of α , respectively, a within nano-spherical domains of Ag, Au, and Pt based on Mindlin’s second strain gradient theory. It is observed that the phenomenon of relaxation for larger spheres (α larger) has little or no effect on the stress field distribution as the center of the sphere is approached, whereas steep variations in the stresses near its surface occur and attain notable values just beneath the surface. Through Gurtin–Murdoch surface elasticity, however, all the non-zero components of the induced stress field in nanosphere are constant and τrr = τθ θ = τφφ . The variation of these stress components in eV/A˚ 3 versus normalized a

sphere radius are plotted in Fig. 8(a)–(c) for Ag, Au, and Pt, respectively. As it is seen, the surface effect is more a0 pronounced for smaller spheres; in all the considered cases the variation of stresses becomes sharper as a/a0 → 0. Likewise, r in the case of an embedded spherical nanocavity, the variation of τ rr and τθ θ = τφφ in eV/A˚ 3 versus for different values of a cavity size are displayed in Figs. 9 and 10 for Ag, Au, and Pt based on Mindlin’s second strain gradient theory and Gurtin– Murdoch surface elasticity, respectively. In Mindlin’s solution, generally, the variation of the stress components near the boundary of the cavity have oscillatory nature, as it was the case in the case of the displacement field. The stresses attain

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

695

Fig. 10. Variation of the stress field components τ rr and τθ θ = τφφ in eV/A˚ 3 versus r/a in the (a) Ag, (b) Au, and (c) Pt domain solid with spherical cavity a . of radius a according to Gurtin–Murdoch surface elasticity for different ratios of α = a0

notably large values at the cavity’s surface, just inside the matrix, and decay rapidly with distance from the cavity. A similar phenomenon is observed within Gurtin–Murdoch surface elasticity, except the solutions are not oscillatory. 8. Conclusion Mindlin’s second strain gradient theory has been formulated in an arbitrary orthogonal curvilinear coordinate system. Using Eringen [12] mathematical tools for transformation from Cartesian coordinates to any arbitrary curvilinear coordinates, the equilibrium equations, generalized stress-strain constitutive relations, components of the generalized strain tensors, and three different types of traction boundary conditions in any orthogonal curvilinear coordinate system are derived. In addition, in order to give a highly-pragmatic example in the field of nanomechanics, Mindlin’s second strain gradient theory is represented in the spherical coordinate system and surface relaxation associated with Ag, Au, and Pt nanospheres as well as nanocavities buried in Ag, Au, and Pt is examined in the framework of Mindlin’s second strain gradient theory. The results are compared with those obtained using Gurtin–Murdoch surface elasticity. For further verification, Ag solid nanosphere has

696

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

also been simulated using molecular dynamics with various EAMs. It is observed that the phenomenon of relaxation captured within Mindlin’s second strain gradient theory is remarkably affected by the size of the nanospherical domain while that of Gurtin-Murdoch surface elasticity is nearly size independent. Appendix A. The expressions τ¯ i j , i, j = r, θ , φ in spherical coordinate system

1 1 1 τ¯rr = τrr − τrrr,r + τrφ r,φ + τrθ r,θ + 2τrrr − τφφ r − τrφφ − τθ θ r − τrθ θ r r sin φ r  1 + τrφ r cot φ + 2 2τrrrr − 4τrrθ θ + 4τrrφ r cot φ − 4τrrφφ − 6τrθ θ r − 2τrθ φθ cot φ r

− 7τrφφ r + τrφφ r cot2 φ − τrφφ r csc2 φ − 2τrφφφ cot φ + 2τθ θ θ θ − 2τθ θ φ r cot φ + 2τθ θ φφ + 2τθ φφθ − 2τφφφ r cot φ + 2τφφφφ + 4τrrφ r,φ − τrθ θ r,φ cot φ − 2τrθ φθ ,φ + 2τrφφ r,φ cot φ − 2τrφφφ ,φ − 2τθ θ φ r,φ − 2τφφφ r,φ + τrφφ r,φφ + 4τrrθ r,θ cscφ − 2τrθ θ θ ,θ cscφ + 2τrθ φ r,θ cot φ cscφ − 2τrθ φφ ,θ cscφ − 2τθ θ θ r,θ cscφ



− 2τθ φφ r,θ cscφ + 2τrθ φ r,θ φ cscφ + τrθ θ r,θ θ csc2 φ +



1 4τrrrr,r − 2τrrθ θ ,r r

+ 2τrrφ r,r cot φ − 2τrrφφ ,r − 3τrθ θ r,r − 3τrφφ r,r + 2τrrφ r,rφ + 2τrrθ r,rθ cscφ



τrrrr,rr , 1 1 1 τ¯rθ = τrθ − τrrθ ,r + τrφθ ,φ + + 2τrrθ + τrθ θ − τφφθ − τθ θ θ τ r r sin φ rθ θ ,θ r

   1 + τrφθ + τrθ φ cot φ + 2 2τrrrθ + 4τrrθ r + 4τrrθ φ cot φ + 4τrrφθ cot φ − 6τrθ θ θ +

r

τrθ θ θ cot2 φ + 4τrθ φ r cot φ − 3τrθ φφ + τrθ φφ cot2 φ + τrθ φφ csc2 φ − 7τrφφθ − 2τθ θ θ r − 2τθ θ θ φ cot φ − 2τθ θ φθ cot φ − 2τθ φφ r − 2τθ φφφ cot φ − 2τφφφθ cot φ + 4τrrφθ ,φ − τrθ θ θ ,φ cot φ + 2τrθ φ r,φ + 2τrθ φφ ,φ cot φ + 2τrφφθ ,φ cot φ − 2τθ θ φθ ,φ − 2τφφφθ ,φ + τrφφθ ,φφ + 4τrrθ θ ,θ cscφ + 2τrθ θ r,θ cscφ + 2τrθ θ φ ,θ cot φ cscφ + 2τrθ φθ ,θ cot φ cscφ − 2τθ θ θ θ ,θ cscφ − 2τθ φφθ ,θ cscφ + 2τrθ φθ ,θ φ cscφ

1 + τrθ θ θ ,θ θ csc2 φ + 4τrrrθ ,r + 2τrrθ r,r + 2τrrθ φ ,r cot φ + 2τrrφθ ,r cot φ − 3τrθ θ θ ,r r

− 3τrφφθ ,r + 2τrrφθ ,rφ + 2τrrθ θ ,rθ cscφ + τrrrθ ,rr , 1 1 1 τ¯rφ = τrφ − τrrφ ,r + τrφφ ,φ + + 2τrrφ + τrφ r − τφφφ − τθ θ φ τ r r sin φ rθ φ ,θ r

   1 + τrφφ − τrθ θ cot φ + 2 4τrrrφ − 8τrrθ θ cot φ + 8τrrφ r + 8τrrφφ cot φ −

2r

− 4τrθ θ r cot φ − 10τrθ θ φ − 2τrθ θ φ cot2 φ + τrθ φθ − 3τrθ φθ cot2 φ − τrθ φθ csc2 φ + 4τrφφ r cot φ − 15τrφφφ + τrφφφ cot2 φ − τrφφφ csc2 φ + 4τθ θ θ θ cot φ − 4τθ θ φ r − 4τθ θ φφ cot φ + 4τθ φφθ cot φ − 4τφφφ r − 4τφφφφ cot φ + 8τrrφφ ,φ − 2τrθ θ φ ,φ cot φ − 4τrθ φθ ,φ cot φ + 4τrφφ r,φ + 4τrφφφ ,φ cot φ − 4τθ θ φφ ,φ − 4τφφφφ ,φ + 2τrφφφ ,φφ + 8τrrθ φ ,θ cscφ − 4τrθ θ θ ,θ cot φ cscφ + 4τrθ φ r,θ cscφ + 4τrθ φφ ,θ cot φ cscφ



− 4τθ θ θ φ ,θ cscφ − 4τθ φφφ ,θ cscφ + 4τrθ φφ ,θ φ cscφ + 2τrθ θ φ ,θ θ csc2 φ +



1 4τrrrφ ,r r

− 2τrrθ θ ,r cot φ + 2τrrφ r,r + 2τrrφφ ,r cot φ − 3τrθ θ φ ,r − 3τrφφφ ,r + 2τrrφφ ,rφ



+ 2τrrθ φ ,rθ cscφ + τrrrφ ,rr ,



1 1 1 τ¯θ r = τθ r − τθ rr,r + τθ φ r,φ + τθ θ r,θ + 3τθ rr − τθ φφ − τθ θ θ + 2τθ φ r cot φ r r sin φ r 1 + 2 6τrrθ r − 6τrθ θ θ θ + 12τrθ φ r cot φ − 6τrθ φφ − (2 + cot2 φ )τθ θ θ r − 4τθ θ φθ cot φ r

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

− 2τθ φφ r + 5τθ φφ r cot2 φ − 3τθ φφ r csc2 φ − 4τθ φφφ cot φ + 6τrθ φ r,φ − τθ θ θ r,φ cot φ − 2τθ θ φθ ,φ + 4τθ φφ r,φ cot φ − 2τθ φφφ ,φ + τθ φφ r,φφ + 6τrθ θ r,θ cscφ − 2τθ θ θ θ ,θ cscφ + 4τθ θ φ r,θ cot φ cscφ − 2τθ θ φφ ,θ cscφ + 2τθ θ φ r,θ φ cscφ + τθ θ θ r,θ θ csc2 φ +





1 6τrrθ r,r − 2τrθ θ θ ,r + 4τrθ φ r,r cot φ − 2τrθ φφ ,r − τθ θ θ r,r − τθ φφ r,r + 2τrθ φ r,rφ r



+ 2τrθ θ r,rθ cscφ + τrrθ r,rr ,



  1 1 1 τ¯θ θ = τθ θ − τθ rθ ,r + τθ φθ ,φ + τθ θ θ ,θ + 3τrθ θ + τθ θ r + 2τφθ θ + τθ θ φ cot φ r r sin φ r 1 + 2 6τrrθ θ + 6τrθ θ r + 6τrθ θ φ cot φ + 12τrθ φθ cot φ − 2τθ θ θ θ csc2 φ − 7/2τθ θ φφ r

+ 3/2τθ θ φφ csc2 φ − 4τθ φφθ + 6τrθ φθ ,φ + 2τθ θ φ r,φ + 1/2 cot φ (12τθ θ φ r + 5τθ θ φφ cot φ + 4τθ φφθ cot φ − 2τθ θ θ θ ,φ + 4τθ θ φφ ,φ + 8τθ φφθ ,φ ) + τθ φφθ ,φφ + 6τrθ θ θ ,θ cscφ + 2τθ θ θ r,θ cscφ + 2τθ θ θ φ ,θ cot φ cscφ + 4τθ θ φθ ,θ cot φ cscφ + 2τθ θ φθ ,θ φ cscφ





1 + τθ θ θ θ ,θ θ csc2 φ + 6τrrθ θ ,r + 2τrθ θ r,r + 2τrθ θ φ ,r cot φ + 4τrθ φθ ,r cot φ − τθ θ θ θ ,r r



− τθ φφθ ,r + 2τrθ φθ ,rφ + 2τrθ θ θ ,rθ cscφ + τrrθ θ ,rr ,

 1 1 1 τ¯θ φ = τθ φ − τθ rφ ,r + τθ φφ ,φ + + 3τrθ φ + τθ φ r + τφθ φ + τθ φφ τ r r sin φ θ θ φ ,θ r  1  − τθ θ θ cot φ + 2 6τrrθ φ − 6τrθ θ θ cot φ + 6τrθ φ r + 12τrθ φφ cot φ − 2τθ θ θ r cot φ r

− τθ θ θ φ − 2τθ θ θ φ cot2 φ − τθ θ φθ − 5τθ θ φθ cot2 φ + τθ θ φθ csc2 φ + 4τθ φφ r cot φ − 2τθ φφφ + 5τθ φφφ cot2 φ − 3τθ φφφ csc2 φ + 6τrθ φφ ,φ − τθ θ θ φ ,φ cot φ − 2τθ θ φθ ,φ + 2τθ φφ r,φ + 4τθ φφφ ,φ cot φ + τθ φφφ ,φφ + 6τrθ θ φ ,θ cscφ − 2τθ θ θ θ ,θ cot φ cscφ + 2τθ θ φ r,θ cscφ + 4τθ θ φφ ,θ cot φ cscφ + 2τθ θ φφ ,θ φ cscφ + τθ θ θ φ ,θ θ csc2 φ +





1 6τrrθ φ ,r − 2τrθ θ θ ,r cot φ + 2τrθ φ r,r + 4τrθ φφ ,r cot φ − τθ θ θ φ ,r − τθ φφφ ,r r



+ 2τrθ φφ ,rφ + 2τrθ θ φ ,rθ cscφ + τrrθ φ ,rr ,

1 1 1 τ¯φ r = τφ r − τφ rr,r + τφφ r,φ + + 2τφ rr + τrφ r − τφφφ − τφθ θ τ r r sin φ φθ r,θ r

   1 + τφφ r − τθ θ r cot φ + 2 6τrrφ r − 6τrθ θ r cot φ − 6τrθ φθ + 6τrφφ r cot φ r

− 6τrφφφ + 2τθ θ θ θ cot φ − 3τθ θ φ r − 5τθ θ φ r cot2 φ + 2τθ θ φ r csc2 φ + 2τθ θ φφ cot φ − 2τθ φφθ cot φ − 3τφφφ r + τφφφ r cot2 φ − τφφφ r csc2 φ − 2τφφφφ cot φ + 6τrφφ r,φ − 3τθ θ φ r,φ cot φ − 2τθ φφθ ,φ + 2τφφφ r,φ cot φ − 2τφφφφ ,φ + τφφφ r,φφ + 6τrθ φ r,θ cscφ − 2τθ θ θ r,θ cot φ cscφ − 2τθ θ φθ ,θ cscφ + 2τθ φφ r,θ cot φ cscφ − 2τθ φφφ ,θ cscφ



+ 2τθ φφ r,θ φ cscφ + τθ θ φ r,θ θ csc2 φ +



1 6τrrφ r,r − 2τrθ θ r,r cot φ − 2τrθ φθ ,r r

+ 2τrφφ r,r cot φ − 2τrφφφ ,r − τθ θ φ r,r − τφφφ r,r + 2τrφφ r,rφ + 2τrθ φ r,rθ cscφ



τrrφ r,rr ,  1 1 1 = τφθ − τφ rθ ,r + τφφθ ,φ + + 3τrφθ + τφθ r + τφθ φ + τφφθ τ r r sin φ φθ θ ,θ r  1  − τθ θ θ cot φ + 2 6τrrφθ − 6τrθ θ θ cot φ + 6τrθ φ r + 6τrθ φφ cot φ + 6τrφφθ cot φ +

τ¯φθ

r

− 2τθ θ θ r cot φ − 2τθ θ θ φ cot2 φ − 3τθ θ φθ − 6τθ θ φθ cot2 φ + 2τθ θ φθ csc2 φ + 4τθ φφ r cot φ − 2τθ φφφ + 2τθ φφφ cot2 φ − 3τφφφθ + 6τrφφθ ,φ − 3τθ θ φθ ,φ cot φ + 2τθ φφ r,φ

697

698

F. Ojaghnezhad and H.M. Shodja / Applied Mathematical Modelling 76 (2019) 669–698

+ 2τθ φφφ ,φ cot φ + 2τφφφθ ,φ cot φ + τφφφθ ,φφ + 6τrθ φθ ,θ cscφ − 2τθ θ θ θ ,θ cot φ cscφ + 2τθ θ φ r,θ cscφ + 2τθ θ φφ ,θ cot φ cscφ + 2τθ φφθ ,θ cot φ cscφ + 2τθ φφθ ,θ φ cscφ





1 + τθ θ φθ ,θ θ csc2 φ + 6τrrφθ ,r − 2τrθ θ θ ,r cot φ + 2τrθ φ r,r + 2τrθ φφ ,r cot φ r



+ 2τrφφθ ,r cot φ − τθ θ φθ ,r − τφφφθ ,r + 2τrφφθ ,rφ + 2τrθ φθ ,rθ cscφ + τrrφθ ,rr ,

 1 1 1 τ¯φφ = τφφ − τφ rφ ,r + τφφφ ,φ + + 3τrφφ + τφφ r + τφφφ − τθ θ φ τ r r sin φ φθ φ ,θ r  1  − τφθ θ cot φ + 2 6τrrφφ − 6τrθ θ φ cot φ − 6τrθ φθ cot φ + 6τrφφ r + 6τrφφφ cot φ r

+ 2τθ θ θ θ cot

φ − 4τθ θ φ r cot φ − 2τθ θ φφ − 6τθ θ φφ cot2 φ + 2τθ θ φφ csc2 φ − τθ φφθ 3τθ φφθ cot2 φ + τθ φφθ csc2 φ + 2τφφφ r cot φ − 3τφφφφ + τφφφφ cot2 φ − τφφφφ csc2 φ 6τrφφφ ,φ − 3τθ θ φφ ,φ cot φ − 2τθ φφθ ,φ cot φ + 2τφφφ r,φ + 2τφφφφ ,φ cot φ + τφφφφ ,φφ 6τrθ φφ ,θ cscφ − 2τθ θ θ φ ,θ cot φ cscφ − 2τθ θ φθ ,θ cot φ cscφ + 2τθ φφ r,θ cscφ

1 2τθ φφφ ,θ cot φ cscφ + 2τθ φφφ ,θ φ cscφ + csc2 φτθ θ φφ ,θ θ + 6τrrφφ ,r 2

− + + +

r − 2τrθ θ φ ,r cot φ − 2τrθ φθ ,r cot φ + 2τrφφ r,r + 2τrφφφ ,r cot φ − τθ θ φφ ,r − τφφφφ ,r



+ 2τrφφφ ,rφ + 2τrθ φφ ,rθ cscφ + τrrφφ ,rr . References [1] J. Guo, Y. Zhao, The size-dependent elastic properties of nanofilms with surface effects, J. App. Phys. 98 (2005) 074306. [2] F. Ojaghnezhad, H.M. Shodja, Surface elasticity revisited in the context of second strain gradient theory, Mech. Mat. 93 (2016) 220–237. [3] F. Ojaghnezhad, H.M. Shodja, A combined first principles and analytical determination of the modulus of cohesion, surface energy, and the additional constants in the second strain gradient elasticity, Int. J. Solids Struct. 50 (2013) 3967–3974. [4] M. Lazar, Dislocations and Cracks in Generalized Continua. Encyclopedia of Continuum Mechanics, Springer-Verleg GmbH, Germany, 2018. [5] M. Lazar, E. Agiasofitu, Fundamentals in generalized elasticity and dislocation theory of quasicrystals: Green tensor, dislocation key-formulas and dislocation loops, Phil. Mag. 94 (35) (2014) 4080–4101. [6] R.A. Toupin, Theories of elasticity with couple-stress, Arch. Ration. Mech. Anal. 17 (1964) 85–112. [7] R.A. Toupin, D.C. Gazis, Surface effects and initial stress in continuum and lattice models of elastic crystals, in: Wallis (Ed.), Proceedings of the International Conference on Lattice Dynamics, Pergamon Press, Copenhagen, 1964, pp. 597–602. August 1963. [8] R.D. Mindlin, Second gradient of strain and surface-tension in linear elasticity, Int. J. Solids Struct. 1 (1965) 417–438. [9] H.M. Shodja, F. Ahmadpoor, A. Tehranchi, Calculation of the additional constants for fee materials in second strain gradient elasticity: Behavior of a nano-size bernoulli-euler beam with surface effects, J. App. Mech. 79 (2012) 0211008. [10] Y. Zhang, Y. Zhao, Measuring the nonlocal effects of a micro/nanobeam by the shifts of resonant frequencies, Int. J. Solids Struct. 102–103 (2016) 259–266. [11] M.E. Gurtin, A.I. Murdoch, A continuum theory of elastic material surfaces, Arch. Rat. Mech. Anal. 57 (1975). [12] A.C. Eringen, Mechanics of Continua, Robert E. Krieger publishing Company, New York, 1980. [13] X. Ji, A.Q. Li, S.J. Zhou., The strain gradient elasticity theory in orthogonal curvilinear coordinates and its applications, J. Mech. (2016) 1–13. [14] S. Zhou, A. Li, B. Wang, A reformulation of constitutive relations in the strain gradient elasticity theory for isotropic materials, Int. J. Solids Struct. 80 (2016) 28–37. [15] J. Grandidier, D.M. Callahan, J.N. Munday, H.A. Atwater, Light absorption enhancement in thin-film solar cells using whispering gallery modes in dielectric nanospheres, Adv. Mater. 23 (2011) 1272–1276. [16] C.E. Roman-Velazquez, C. Noguez, J.Z. Zhang, Theoretical study of surface plasmon resonances in hollow gold-silver double-shell nanostructures, J. Phys. Chem. 113 (2009) 4068–4074. [17] S.W. Kim, M. Kim, W.Y. Lee, T. Hyeon, Fabrication of hollow palladium spheres and their successful application to the recyclable heterogeneous catalyst for suzuki coupling reactions, J. Am. Chem. Soc. 124 (26) (2002) 7642–7643. [18] H.P. Liang, L.J. Wan, C.L. Bai, L. Jiang, Gold hollow nanospheres: Tunable surface plasmon resonance controlled by interior-cavity sizes, J. Phys. Chem. B. 109 (16) (2005) 7795–7800. [19] W. Lu, G.D. Zhang, R. Zhang, L.G. Flores, Q. Huang, J.G. Gelovani, C. Li, Tumor site-specific silencing of NF-kappab p65 by targeted hollow gold nanosphere-mediated photothermal transfection, Cancer Res. 70 (8) (2010) 3177–3188. [20] L.R. Hirsch, R.J. Stafford, J.A. Bankson, S.R. Sershen, B. Rivera, R.E. Price, J.D. Hazle, N.J. Halas, J.L. West, Nanoshell-mediated near-infrared thermal therapy of tumors under magnetic resonance guidance, Proc. Natl. Acad. Sci. 100 (23) (2003) 13549–13554. [21] M.E. Gurtin, A.I. Murdoch, Surface stress in solids, Int. J. Solids Struct. 14 (1978) 431–440. [22] H.M. Shodja, H. Moosavian, F. Ojaghnezhad, Toupin-Mindlin first strain gradient theory revisited for cubic crystals of hexoctahedral class: analytical expression of the material parameters in terms of the atomic force constants and evaluation via ab initio DFT, Mech. Mat. 123 (2018) 19–29. [23] S.M. Foiles, M.I. Baskes, M.S. Daw, Embedded-atom-method functions for the FCC metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys, Phys. Rev. B. 33 (1986) 7983. [24] G.J. Ackland, G. Tichy, V. Vitek, M.W. Finnis, Simple N-body potentials for the noble metals and nickel, Phil. Mag. A. 56 (6) (1986) 735–756. [25] P.L. Williams, Y. Mishin, J.C. Hamilton, An embedded-atom potential for the Cu-Ag system, Model. Simul. Mater. Sci. Eng. 14 (2006) 817–833. [26] X.W. Zhou, R.A. Johnson, H.N.G. Wadley, Misfit-energy-increasing dislocations in vapor-deposited cofe/nife multilayers, Phys. Rev. B 69 (2004) 144113.