Axisymmetric loading on nanoscale multilayered media

Abstract

Multilayered nanoscale structures are used in several applications. Because the effect of surface energy becomes nontrivial at such a small scale, a modified continuum theory is required to accurately predict their mechanical behaviors. A Gurtin–Murdoch continuum model of surface elasticity is implemented to establish a computational scheme for investigating an elastic multilayered system under axisymmetric loads with the incorporation of surface/interface energy. Each layer stiffness matrix is derived based on the general solutions of stresses and displacements obtained in the form of the Hankel integral transform. Numerical solutions to the global equation, which are formulated based on the continuity conditions of tractions and displacements across interfaces between layers, yield the displacements at each layer interface and on the top surface of the multilayered medium. The numerical solutions indicate that the elastic responses of multilayered structures are affected significantly by the surface material properties of both the top surface and interfaces, and that they become size dependent. In addition, the indentation problem of a multilayered nanoscale elastic medium under a rigid frictionless cylindrical punch is investigated to demonstrate the application of the proposed solution scheme.

Graphical abstract

Keywords

functionally graded layer / Gurtin–Murdoch surface elasticity / multilayered medium / size dependency / stiffness matrix

1 Introduction

In recent decades, nanotechnology has been applied extensively in diverse fields, such as engineering, manufacturing, medicine, and energy production. Owing to their superior mechanical, electronic, and optical properties, nanoscale structures have been developed for applications in electronic conductors, magnetic storage media, and hard surface coatings. The mechanical behavior of nanoscale structures (i.e., structures that measure between 1 and 100 nm) become size dependent and cannot be precisely described by the conventional continuum theory [1,2]. A comprehensive and profound understanding of such complex behaviors is crucial for the development of nanoscale devices and structures.

Although experimental studies have helped determine the actual physical behaviors of structures at the nanoscale, they generally require high-precision testing instruments and procedures. Atomistic and molecular simulations are effective techniques for accurately characterizing the mechanical responses of nanoscale structures. However, they are typically computationally prohibitive in practical applications as they require significant computational effort for modeling the numerous atoms and molecules associated with the systems. Hence, a modified continuum-based approach should be devised to account for the effect of nanoscale structures in obtaining first approximations for the examination of nanoscale structures.

Regarding nanoscale structures, the contribution of surface/interface free energy effects from atoms near the surface/interface is nontrivial, and their mechanical behaviors are size dependent owing to their high surface-to-volume ratio. To account for these effects, Gurtin and Ian Murdoch (G–M) formulated a surface elasticity theory based on the thermodynamic principles of elastic material surfaces [3,4]. The model was validated using atomistic models of bars and plates at the nanoscale, under both uniaxial loading and pure bending [5]. In their study [5], the surface was regarded as an extremely thin membrane that completely adhered to the bulk material, and its elastic modulus was obtained from atomistic simulations [5,6]. In addition, the effect of the surface energy was considerable for extremely soft materials, such as polymer gels. Over the past two decades, the G–M model [3,4] has been successfully employed to examine the size-dependent mechanical behaviors of various boundary value problems, including those involving ultra-thin elastic films [7], layered elastic media [8], thin plates and nanoplates [9,10], nanostructure optimization [11], nanocracks [12,13], nanoinhomogeneities [14], and nanoindentations [15–18].

Various continuum-based concepts have been proposed to describe microstructural effects and size-dependent behaviors such as nonlocal elasticity [19,20], couple stress elasticity [21,22], and strain gradient elasticity [23], which are typically overlooked in macroscopic continuum theory. Wongviboonsin et al. [24] investigated the effect of material microstructure on the elastic responses of a layered medium based on coupled stress elasticity. Efficient numerical solutions based on neural network methods have been proposed recently for solving boundary value problems (e.g., [25,26]); these solutions do not require discretization and serve as a natural framework for inverse analysis and optimization. To consider size-dependent behaviors, Nguyen-Thanh et al. [27] presented machine learning-based solutions for solving continuum mechanics problems by considering the strain gradient elasticity. Notably, machine learning- and numerical technique-based solutions, such as boundary element method (BEM) and finite element method (FEM), are not the focus of this study.

The mechanical behaviors of multilayered media must be investigated to promote the further development of various applications, such as electrical circuits in electrical engineering, micro- and nanoelectromechanical systems (MEMS/NEMS), and composite laminate materials used in aerospace engineering and construction [28]. Natural soil has been investigated previously for solving various geomechanics problems, since the profiles of natural soil are multilayered (e.g., [29–33]). Meanwhile, multilayered structures are present in various products adopting nanoscale structures, such as semiconductor devices, magnetic/optical data storage devices, dielectric and insulating materials, and hard surface coatings [34,35]. However, nanoscale effects on multilayered structures have not been comprehensively investigated. Zhu et al. [34] studied the mechanical behaviors of nanoscale Cu/Ni multilayers based on nanoindentation, X-ray diffraction, and transmission electron microscopy, which shows that a profound understanding of the elastic fields of multilayered systems is essential for the development of nanoscale devices and systems. To the best of our knowledge, the examination of three-dimensional multilayered media under the effects of surface/interface energy, including in cases involving axisymmetric deformations, has not yet been reported.

In the current study, a semi-analytical solution scheme for investigating the effects of surface energy on nanoscale multilayered structures under axisymmetric deformations is developed based on the G–M model and the exact stiffness matrix method [36]. Surface energy effects are considered for both the top surface and at each interface between the layers. The relevant stiffness matrix is obtained from general solutions expressed in the form of the Hankel integral transform [37]. The proposed solution method is validated via a comparison with the available solutions. Several numerical results are illustrated for bi-material multilayered and functionally graded (FG) layers to demonstrate the effects of surface/interface energy. Additionally, a frictionless indentation problem of a bi-material multilayered system under a flat-ended cylindrical indenter is examined to demonstrate the application of the present solution scheme. The solution scheme utilized in the current study allows one to gain fundamental understanding pertaining to the mechanical characteristics of multilayered structures under the effects of surface/interface energy. Additionally, the scheme can be employed as a benchmark for the further development of numerical schemes such as the BEM and FEM.

2 Methodology

2.1 Problem description

A three-dimensional elastic medium overlying a rigid base, comprising N layers with various properties and thicknesses, under a cylindrical coordinate system and governed by the G–M surface elasticity [3,4], is shown in Fig.1. Each layer of the multilayered medium is separated into two sections, i.e., the bulk medium and a zero-thickness surface membrane attached completely to the bulk. The bulk of the nth layer (n = 1,2,…,N) is an isotropic elastic solid with Lamé constants λ⁽ⁿ⁾ and μ⁽ⁿ⁾, whereas the surface n (n = 1,2,…,N) possesses surface Lamé constants {\lambda }_n^{\mathrm{s}} and {\mu }_n^{\mathrm{s}}, and a residual surface stress {\tau }_n^{\mathrm{s}} (assumed to be constant). The superscript (n) denotes the layer number, whereas the subscript n denotes the surface number. The nth surface with n > 1 can be considered an interface that represents the merged surface from the bottom surface of the upper layer to the top surface of the lower layer. In addition, the multilayered medium is imposed with only an axisymmetric vertical load applied on its top surface, denoted by p(r), with no remote loading or body force.

The axisymmetric field equations and the corresponding general solutions of the bulk of each layer can be obtained based on the classical linear elasticity. For the nth surface (n = 1,2,…,N), the field equations under the axisymmetric condition are expressed by the generalized Young–Laplace equations [38], the stress–strain equations, and the strain–displacement equations of the surface/interface, as follows [3,4,39].

(1)\begin{array}{rl}& \frac{\mathrm{\partial }{\sigma }_{r\mathit{\text{z}}n}^{\mathrm{s}}}{\mathrm{\partial }r}+\frac{{\sigma }_{r\mathit{\text{z}}n}^{\mathrm{s}}}{r}+t_{\mathit{\text{z}}n}^{\mathrm{t}}+t_{\mathit{\text{z}}n}^{\mathrm{b}}=0,\\ & \frac{\mathrm{\partial }{\sigma }_{rrn}^{\mathrm{s}}}{\mathrm{\partial }r}+\frac{{\sigma }_{rrn}^{\mathrm{s}}-{\sigma }_{\theta \theta n}^{\mathrm{s}}}{r}+t_{rn}^{\mathrm{t}}+t_{rn}^{\mathrm{b}}=0,\end{array}

\begin{array}{}\text{(2a)}& {\sigma }_{rrn}^{\mathrm{s}}=\left({\lambda }_n^{\mathrm{s}}+2{\mu }_n^{\mathrm{s}}\right){\epsilon }_{rrn}^{\mathrm{s}}+\left({\lambda }_n^{\mathrm{s}}+{\tau }_n^{\mathrm{s}}\right){\epsilon }_{\theta \theta n}^{\mathrm{s}}+{\tau }_n^{\mathrm{s}},\end{array}

\begin{array}{}\text{(2b)}& {\sigma }_{\theta \theta n}^{\mathrm{s}}=\left({\lambda }_n^{\mathrm{s}}+2{\mu }_n^{\mathrm{s}}\right){\epsilon }_{\theta \theta n}^{\mathrm{s}}+\left({\lambda }_n^{\mathrm{s}}+{\tau }_n^{\mathrm{s}}\right){\epsilon }_{rrn}^{\mathrm{s}}+{\tau }_n^{\mathrm{s}},\end{array}

\begin{array}{}\text{(2c)}& {\sigma }_{r\mathit{\text{z}}n}^{\mathrm{s}}={\tau }_n^{\mathrm{s}}\frac{\mathrm{d}{u}_{\mathit{\text{z}}n}^{\mathrm{s}}}{\mathrm{d}r},\end{array}

(3){\epsilon }_{rrn}^{\mathrm{s}}=\frac{\mathrm{d}{u}_{rn}^{\mathrm{s}}}{\mathrm{d}r},{\epsilon }_{\theta \theta n}^{\mathrm{s}}=\frac{u_{rn}^{\mathrm{s}}}{r}.

In the equations above, the superscript “s” represents the parameters associated with the surface; t_{\mathit{\text{z}}n}^{\mathrm{t}} and t_{rn}^{\mathrm{t}} denote the vertical and radial tractions imposed on the top side of the surface n, respectively; t_{\mathit{\text{z}}n}^{\mathrm{b}} and t_{rn}^{\mathrm{b}} represent the vertical and radial tractions imposed on the bottom side of the surface n, respectively; u_{\mathit{\text{z}}n}^{\mathrm{s}} and u_{rn}^{\mathrm{s}} denote the vertical and radial displacements of the surface n, respectively. Additionally, {\sigma }_{rrn}^{\mathrm{s}} and {\sigma }_{\theta \theta n}^{\mathrm{s}} represent the in-plane components of surface stresses of the surface n, whereas {\sigma }_{r\mathit{\text{z}}n}^{\mathrm{s}} denotes the out-of-plane component of surface stresses of the surface n.

The boundary conditions imposed at the top surface and rigid base are as follows.

\begin{array}{}\text{(4a)}& \left\{{\sigma }_{\mathit{\text{zz}}}^{\left(1\right)}+T_{\mathit{\text{z}}}^{\mathrm{s}}\right\}|\genfrac{}{}{0pt}{}{}{\mathit{\text{z}}={\mathit{\text{z}}}_1}=-p\left(r\right),\end{array}

\begin{array}{}\text{(4b)}& \left\{{\sigma }_{\mathit{\text{z}}r}^{\left(1\right)}+T_r^{\mathrm{s}}\right\}|\genfrac{}{}{0pt}{}{}{\mathit{\text{z}}={\mathit{\text{z}}}_1}=0,\end{array}

\begin{array}{}\text{(4c)}& u_i^{\left(N\right)}|\genfrac{}{}{0pt}{}{}{\mathit{\text{z}}={\mathit{\text{z}}}_{N+1}}=0,i=r,\mathit{\text{z}},\end{array}

where

\begin{array}{}\text{(5a)}& T_{\mathit{\text{z}}}^{\mathrm{s}}|\genfrac{}{}{0pt}{}{}{\mathit{\text{z}}={\mathit{\text{z}}}_1}={\tau }_1^{\mathrm{s}}{\left(\frac{{\mathrm{d}}^2{u}_{\mathit{\text{z}}}^{\mathrm{s}}}{\mathrm{d}{r}^2}+\frac{1}{r}\frac{\mathrm{d}{u}_{\mathit{\text{z}}}^{\mathrm{s}}}{\mathrm{d}r}\right)}_{\mathit{\text{z}}={\mathit{\text{z}}}_1},\end{array}

\begin{array}{}\text{(5b)}& T_r^{\mathrm{s}}|\genfrac{}{}{0pt}{}{}{\mathit{\text{z}}={\mathit{\text{z}}}_1}={\kappa }_1^{\mathrm{s}}{\left(\frac{{\mathrm{d}}^2{u}_r^{\mathrm{s}}}{\mathrm{d}{r}^2}+\frac{1}{r}\frac{\mathrm{d}{u}_r^{\mathrm{s}}}{\mathrm{d}r}-\frac{u_r^{\mathrm{s}}}{r^2}\right)}_{\mathit{\text{z}}={\mathit{\text{z}}}_1},\end{array}

where {\kappa }_1^{\mathrm{s}}={\lambda }_1^{\mathrm{s}}+2{\mu }_1^{\mathrm{s}}. The variables T_{\mathit{\text{z}}}^{\mathrm{s}} and T_r^{\mathrm{s}}, which do not exist in classical elasticity, represent the contributions of surface energy effects in the tangential and normal directions, respectively. Both T_{\mathit{\text{z}}}^{\mathrm{s}} and T_r^{\mathrm{s}} occur at the top surface and at the interfaces of the multilayered medium under the effect of surface energy. Notably, T_{\mathit{\text{z}}}^{\mathrm{s}} and T_r^{\mathrm{s}} in Eq. (5) are expressed based on the assumption of a constant residual surface stress, i.e., {\displaystyle \frac{\mathrm{\partial }{\tau }_1^{\mathrm{s}}}{\mathrm{\partial }r}}=0. At each interface (the nth surface with n = 2,3,…,N), the traction and displacement continuity conditions can be expressed as follows, respectively.

\begin{array}{}\text{(6a)}& \left\{{\sigma }_{\mathit{\text{z}}i}^{(n-1)}-{\sigma }_{\mathit{\text{z}}i}^{\left(n\right)}-T_i^{\mathrm{s}}\right\}|\genfrac{}{}{0pt}{}{}{\mathit{\text{z}}={\mathit{\text{z}}}_n}=0,i=r,\mathit{\text{z}},\end{array}

\begin{array}{}\text{(6b)}& u_i^{(n-1)}|\genfrac{}{}{0pt}{}{}{\mathit{\text{z}}={\mathit{\text{z}}}_n}=u_i^{\mathrm{s}}|\genfrac{}{}{0pt}{}{}{\mathit{\text{z}}={\mathit{\text{z}}}_n}=u_i^{\left(n\right)}|\genfrac{}{}{0pt}{}{}{\mathit{\text{z}}={\mathit{\text{z}}}_n},i=r,\mathit{\text{z}}.\end{array}

Moreover, the following far-field regularity condition is satisfied in the bulk of each layer.

\begin{array}{}\text{(6c)}& u_i(r,\mathit{\text{z}})\to 0,r\to \mathrm{\infty },i=r,\mathit{\text{z}}.\end{array}

When an axisymmetric radial or vertical loading is exerted at the nth interface, the right-hand side of Eq. (6a) with i = r (for the radial load) or z (for the vertical load) will not be zero; in fact, it will be expressed in terms of the relevant loading.

2.2 Solution of boundary value problem

To obtain a solution to the boundary value problem shown in Fig.1, the boundary conditions expressed in Eqs. (4a)–(4c) are considered along with the continuity conditions presented in Eqs. (6a) and (6b), respectively. Because the problem is rotationally symmetric about the z-axis, the Hankel integral transform method can be conveniently employed [37]. A solution scheme based on the conventional method, which directly solves the boundary value problem using Eqs. (4) and (6), yields a high condition number for the equation system when the value of the Hankel transform parameter increases, owing to both negative and positive exponential terms appearing in the equations [36]. By adopting the stiffness matrix approach [36], the equations will not contain unstable exponential functions, which results in accurate numerical computations. Thus, the stiffness matrix scheme was adopted in this study to investigate the multilayered medium under surface energy effects.

According to Sneddon [37], displacements that satisfy the far-field regularity condition shown in Eq. (6c) at the upper and lower surfaces of the bulk in the nth layer of the multilayered medium shown in Fig.1 are expressed in the Hankel integral transform space as follows.

\begin{array}{}\text{(7a)}& \left\{\begin{array}{c}{\overline{u}}_{\mathit{\text{z}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ {\overline{u}}_r(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ {\overline{u}}_{\mathit{\text{z}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\\ {\overline{u}}_r(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\end{array}\right\}=\left[\begin{array}{c}{\mathit{G}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ ...\\ {\mathit{G}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\end{array}\right]{\mathit{C}}^{\left(n\right)}\left(\overline{\xi }\right),\end{array}

where ξ represents the Hankel transform parameter, and

\begin{array}{}\text{(7b)}& {\mathit{C}}^{\left(n\right)}\left(\overline{\xi }\right)={\left[A^{\left(n\right)}{B}^{\left(n\right)}{C}^{\left(n\right)}{D}^{\left(n\right)}\right]}^{\mathrm{T}},\end{array}

\begin{array}{rl}& {\mathit{G}}^{\left(n\right)}\left(\overline{\xi },\overline{\mathit{\text{z}}}\right)=\overline{\xi }{\mathrm{e}}^{-\overline{\xi }\overline{\mathit{\text{z}}}}\left(\frac{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}{{\overline{\mu }}^{\left(n\right)}}\right)\\ & \cdot \left[\begin{array}{cccc}-\overline{\xi }& -\left({\displaystyle \frac{2{\overline{\mu }}^{\left(n\right)}}{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}}+\overline{\xi }\overline{\mathit{\text{z}}}\right)& -\overline{\xi }{\mathrm{e}}^{2\overline{\xi }\overline{\mathit{\text{z}}}}& \left({\displaystyle \frac{2{\overline{\mu }}^{\left(n\right)}}{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}}-\overline{\xi }\overline{\mathit{\text{z}}}\right){\mathrm{e}}^{2\overline{\xi }\overline{\mathit{\text{z}}}}\\ -\overline{\xi }& \left(1-\overline{\xi }\overline{\mathit{\text{z}}}\right)& \overline{\xi }{\mathrm{e}}^{2\overline{\xi }\overline{\mathit{\text{z}}}}& \left(1+\overline{\xi }\overline{\mathit{\text{z}}}\right){\mathrm{e}}^{2\overline{\xi }\overline{\mathit{\text{z}}}}\end{array}\right],\end{array}

where A^{(n)}, B^{(n)}, C^{(n)}, and D^{(n)} are arbitrary functions corresponding to the bulk of layer n, which can be obtained by imposing boundary and continuity conditions. Additionally,

\begin{array}{}\text{(7d)}& {\overline{\mu }}^{\left(n\right)}=\frac{{\mu }^{\left(n\right)}}{{\mu }^{\left(1\right)}},{\overline{\lambda }}^{\left(n\right)}=\frac{{\lambda }^{\left(n\right)}}{{\mu }^{\left(1\right)}},\end{array}

\begin{array}{}\text{(7e)}& \overline{\mathit{\text{z}}}=\mathit{\text{z}}/\mathrm{\Lambda },\overline{\xi }=\xi \mathrm{\Lambda },\mathrm{\Lambda }=\frac{{\kappa }_1^s\left({\lambda }^{\left(1\right)}+2{\mu }^{\left(1\right)}\right)}{2{\mu }^{\left(1\right)}\left({\lambda }^{\left(1\right)}+{\mu }^{\left(1\right)}\right)},\end{array}

where Λ represents the reference length-scale parameter defined by the elastic properties of the top layer and the material constants of the top surface.

Similarly, the Hankel transforms of the tractions at the overlying and underlying surfaces of the bulk in the nth layer are expressed as

\begin{array}{}\text{(8a)}& \left\{\begin{array}{c}-{\overline{\sigma }}_{\mathit{\text{zz}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ -{\overline{\sigma }}_{\mathit{\text{z}}r}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ {\overline{\sigma }}_{\mathit{\text{zz}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\\ {\overline{\sigma }}_{\mathit{\text{z}}r}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\end{array}\right\}=\left[\begin{array}{c}-{\mathit{H}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ ...\\ {\mathit{H}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\end{array}\right]{\mathit{C}}^{\left(n\right)}\left(\overline{\xi }\right),\end{array}

where

\begin{array}{rl}& {\mathit{H}}^{\left(n\right)}\left(\overline{\xi },\overline{\mathit{\text{z}}}\right)=2{\overline{\xi }}^2{\mathrm{e}}^{-\overline{\xi }\overline{\mathit{\text{z}}}}\left({\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}\right)\\ & \cdot \left[\begin{array}{cccc}\overline{\xi }& \left({\displaystyle \frac{{\overline{\mu }}^{\left(n\right)}}{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}}+\overline{\xi }\overline{\mathit{\text{z}}}\right)& -\overline{\xi }{\mathrm{e}}^{2\overline{\xi }\overline{\mathit{\text{z}}}}& \left({\displaystyle \frac{{\overline{\mu }}^{\left(n\right)}}{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}}-\overline{\xi }\overline{\mathit{\text{z}}}\right){\mathrm{e}}^{2\overline{\xi }\overline{\mathit{\text{z}}}}\\ \overline{\xi }& -\left({\displaystyle \frac{{\overline{\lambda }}^{\left(n\right)}}{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}}-\overline{\xi }\overline{\mathit{\text{z}}}\right)& \overline{\xi }{\mathrm{e}}^{2\overline{\xi }\overline{\mathit{\text{z}}}}& \left({\displaystyle \frac{{\overline{\lambda }}^{\left(n\right)}}{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}}+\overline{\xi }\overline{\mathit{\text{z}}}\right){\mathrm{e}}^{2\overline{\xi }\overline{\mathit{\text{z}}}}\end{array}\right].\end{array}

Based on Eq. (5), T_{\mathit{\text{z}}}^{\mathrm{s}} and T_r^{\mathrm{s}} at the nth interface are expressed in the Hankel transform space as

\begin{array}{}\text{(9a)}& \left\{\begin{array}{c}{\overline{T}}_{\mathit{\text{z}}}^{\mathrm{s}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ {\overline{T}}_r^{\mathrm{s}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\end{array}\right\}={\mathit{S}}_n(\overline{\xi },{\overline{\mathit{\text{z}}}}_n){\mathit{C}}^{\left(n\right)}\left(\overline{\xi }\right),\end{array}

where

\begin{array}{rl}& {\mathit{S}}_n(\overline{\xi },\overline{\mathit{\text{z}}})={\overline{\xi }}^2{\mathrm{e}}^{-\overline{\xi }\overline{\mathit{\text{z}}}}\\ & \cdot \left[\begin{array}{cccc}-\alpha {\overline{\xi }}^2& -\alpha \left({\displaystyle \frac{2\overline{\xi }{\overline{\mu }}^{\left(n\right)}}{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}}+{\overline{\xi }}^2\overline{\mathit{\text{z}}}\right)& -\alpha {\overline{\xi }}^2{\mathrm{e}}^{2\overline{\xi }\overline{\mathit{\text{z}}}}& \alpha \left({\displaystyle \frac{2\overline{\xi }{\overline{\mu }}^{\left(n\right)}}{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}}-{\overline{\xi }}^2\overline{\mathit{\text{z}}}\right){\mathrm{e}}^{\overline{2\xi }\overline{\mathit{\text{z}}}}\\ -\beta {\overline{\xi }}^2& \beta \overline{\xi }\left(1-\overline{\xi }\overline{\mathit{\text{z}}}\right)& \beta {\overline{\xi }}^2{\mathrm{e}}^{\overline{2\xi }\overline{\mathit{\text{z}}}}& \beta \overline{\xi }\left(1+\overline{\xi }\overline{\mathit{\text{z}}}\right){\mathrm{e}}^{\overline{2\xi }\overline{\mathit{\text{z}}}}\end{array}\right],\end{array}

\begin{array}{}\text{(9c)}& \alpha =\left(\frac{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}{{\overline{\mu }}^{\left(n\right)}}\right){\overline{\tau }}_n^{\mathrm{s}},\beta =\left(\frac{{\overline{\lambda }}^{\left(n\right)}+{\overline{\mu }}^{\left(n\right)}}{{\overline{\mu }}^{\left(n\right)}}\right){\overline{\kappa }}_n^{\mathrm{s}},\end{array}

\begin{array}{}\text{(9d)}& {\overline{\tau }}_n^{\mathrm{s}}=\frac{{\tau }_n^{\mathrm{s}}}{{\mu }^{\left(1\right)}\mathrm{\Lambda }},{\overline{\kappa }}_n^{\mathrm{s}}=\frac{{\kappa }_n^{\mathrm{s}}}{{\mu }^{\left(1\right)}\mathrm{\Lambda }},\end{array}

where α and β represent the contributions of the surface/interface energy from the residual surface stress and surface material parameter, respectively.

Based on Eq. (7a), {\mathit{C}}^{\left(n\right)}\left(\overline{\xi }\right) can be written as

(10){\mathit{C}}^{\left(n\right)}\left(\overline{\xi }\right)={\left[\begin{array}{c}{\mathit{G}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ ...\\ {\mathit{G}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\end{array}\right]}^{-1}\left\{\begin{array}{c}{\stackrel{}{u}}_{\mathit{\text{z}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ {\overline{u}}_r(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ {\overline{u}}_{\mathit{\text{z}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\\ {\overline{u}}_r(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\end{array}\right\}.

Substituting the expression above for {\mathit{C}}^{\left(n\right)}\left(\overline{\xi }\right) into Eqs. (8a) and (9a) along with the condition that the tractions T_r^{\mathrm{s}} and T_{\mathit{\text{z}}}^{\mathrm{s}} at the nth interface are imposed at the upper surface of layer n (z = z_n), the relationship between the tractions and displacements at each interface in the nth layer (n = 1,2,…,N) can be written in the following matrix equation.

\begin{array}{}\text{(11a)}& {\mathit{F}}^{\left(n\right)}={\mathit{K}}^{\left(n\right)}{\mathit{U}}^{\left(n\right)},\end{array}

where

\begin{array}{rl}{\mathit{F}}^{\left(n\right)}=& [-\left\{{\overline{\sigma }}_{\mathit{\text{z}}z}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)+{\overline{T}}_{\mathit{\text{z}}}^{\mathrm{s}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\right\}-\left\{{\overline{\sigma }}_{\mathit{\text{z}}r}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)+{\overline{T}}_r^{\mathrm{s}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\right\}\\ & {{\overline{\sigma }}_{\mathit{\text{z}}z}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1}){\overline{\sigma }}_{\mathit{\text{z}}r}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})]}^{\mathrm{T}},\end{array}

\begin{array}{}\text{(11c)}& {\mathit{U}}^{\left(n\right)}={\left[\begin{array}{cccc}{\overline{u}}_{\mathit{\text{z}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)& {\overline{u}}_r(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)& {\overline{u}}_{\mathit{\text{z}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})& {\overline{u}}_r(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\end{array}\right]}^{\mathrm{T}},\end{array}

\begin{array}{}\text{(11d)}& {\mathit{K}}^{\left(n\right)}=\left[\begin{array}{c}-{\mathit{H}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)-{\mathit{S}}_n(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ ...\\ {\mathit{H}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\end{array}\right]{\left[\begin{array}{c}{\mathit{G}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_n)\\ ...\\ {\mathit{G}}^{\left(n\right)}(\overline{\xi },{\overline{\mathit{\text{z}}}}_{n+1})\end{array}\right]}^{-1}.\end{array}

By assembling the stiffness matrix {\mathit{K}}^{\left(n\right)} (n = 1,2,…,N) (Eq. (11a)) using boundary conditions (Eq. (4)), the continuity of the tractions (Eq. (6a)), and the displacements (Eq. (6b)) at all interfaces of the multilayered medium, the global equation system corresponding to the multilayer elastic medium under a vertical applied load and the contribution from the surface energy effects, as shown in Fig.1, can be obtained as follows.

\begin{array}{}\text{(12a)}& {\mathit{F}}_{\mathrm{G}}={\mathit{K}}_{\mathrm{G}}{\mathit{U}}_{\mathrm{G}},\end{array}

where

\begin{array}{}\text{(12b)}& {\mathit{F}}_{\mathrm{G}}={\left[\begin{array}{ccccccc}\stackrel{}{p}\left(\overline{\xi }\right)& 0& 0& 0& \cdots & 0& 0\end{array}\right]}^{\mathrm{T}},\end{array}

\begin{array}{}\text{(12c)}& {\mathit{U}}_{\mathrm{G}}={\left[\begin{array}{ccccccc}{\overline{u}}_{\mathit{\text{z}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_1)& {\overline{u}}_r(\overline{\xi },{\overline{\mathit{\text{z}}}}_1)& {\overline{u}}_{\mathit{\text{z}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_2)& {\overline{u}}_r(\overline{\xi },{\overline{\mathit{\text{z}}}}_2)& \cdots & {\overline{u}}_{\mathit{\text{z}}}(\overline{\xi },{\overline{\mathit{\text{z}}}}_N)& {\overline{u}}_r(\overline{\xi },{\overline{\mathit{\text{z}}}}_N)\end{array}\right]}^{\mathrm{T}},\end{array}

\begin{array}{}\text{(12d)}& {\mathit{K}}_{\mathrm{G}}=\left[\begin{array}{ccccccc}{K}_{11}^{\left(1\right)}& K_{12}^{\left(1\right)}& K_{13}^{\left(1\right)}& K_{14}^{\left(1\right)}& \cdots & 0& 0\\ K_{21}^{\left(1\right)}& K_{22}^{\left(1\right)}& K_{23}^{\left(1\right)}& K_{24}^{\left(1\right)}& \cdots & 0& 0\\ K_{31}^{\left(1\right)}& K_{32}^{\left(1\right)}& K_{33}^{\left(1\right)}+K_{11}^{\left(2\right)}& K_{34}^{\left(1\right)}+K_{12}^{\left(2\right)}& \cdots & 0& 0\\ K_{41}^{\left(1\right)}& K_{42}^{\left(1\right)}& K_{43}^{\left(1\right)}+K_{21}^{\left(2\right)}& K_{44}^{\left(1\right)}+K_{22}^{\left(2\right)}& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & K_{33}^{\left(n-1\right)}+K_{11}^{\left(n\right)}& K_{34}^{\left(n-1\right)}+K_{12}^{\left(n\right)}\\ 0& 0& 0& 0& \cdots & K_{43}^{\left(n-1\right)}+K_{21}^{\left(n\right)}& K_{44}^{\left(n-1\right)}+K_{22}^{\left(n\right)}\end{array}\right],\end{array}

where K_G is a 2N × 2N global stiffness matrix. In addition, the function \stackrel{}{p}\left(\overline{\xi }\right) denotes the Hankel transform of the normalized vertical load p\left(r\right)/{\mu }^{\left(1\right)}.

By solving Eq. (12a), the Hankel transforms of the vertical and radial displacements at the top surface and at all interfaces can be obtained. Substituting the displacement solutions into Eq. (11a) results in the Hankel transforms of tractions at the top and bottom surfaces of the bulk of all the layers. Subsequently, applying an accurate numerical quadrature scheme yields the elastic fields of the multilayered medium. The next section presents a brief outline of the numerical integration scheme employed, followed by its verification based on existing solutions. Subsequently, selected numerical results of displacements and stresses for bi-material systems and FG elastic media are presented.

3 Results and discussion

In this section, numerical solutions based on the proposed scheme are presented to illustrate the effects of surface/interface energy on nanoscale multilayered systems. The proposed solution scheme is first validated using existing benchmark solutions. Subsequently, parametric studies regarding multilayered and FG layer media under surface loadings are performed to illustrate the effects of various parameters, including the surface and size-dependent effects on their elastic fields. Finally, a numerical solution to an indentation problem under the effect of surface energy is presented to demonstrate the application of the proposed solution technique. In the numerical study, the relevant parameters for the layer materials and surface/interface properties are based on those from previous studies [5,40]. Notably, the sensitivity and uncertainty analyses of the input parameters [41,42] are beyond the scope of this study.

3.1 Numerical scheme and verification

The first task of the proposed numerical scheme is to generate a stiffness matrix corresponding to each layer for a specified value of the Hankel transform parameter ξ. Subsequently, a global stiffness equation in the form of Eq. (12a) is established and solved for the global displacement vector shown in Eq. (12c), for the value of ξ. Next, the displacement solution is substituted into Eq. (11) to obtain the Hankel transforms of traction at the top and bottom surfaces of the bulk of all the layers. Finally, the displacements and tractions in the multilayered medium are obtained from the numerical inversion of the Hankel transform. In this study, a semi-infinite integral was computed using an adaptive numerical integration based on the 21-point Gauss–Kronrod rule [43].

For verification, the solutions obtained by Katebi and Selvadurai [44] and Tirapat et al. [8] were compared with the present solutions. Katebi and Selvadurai [44] presented the solutions for an FG layer with a thickness h bonded to an underlying homogeneous half-space under a uniform vertical patch load of radius a applied over the interface. In their model, an exponential variation along the thickness direction of the shear modulus, expressed as \mu \left(\mathit{\text{z}}\right)={\mu }_0{\mathrm{e}}^{\bar{m}\mathit{\text{z}}}, was assumed. In the expression, \bar{m} is a constant and {\mu }_0 is the shear modulus at the top surface of the FG layer; meanwhile, the Poisson’s ratio was set as 0.5. In the current study, the FG layer was modeled by equally discretizing a finite layer into several sublayers with a constant shear modulus specified using the exponential function above at mid-height. The radial profiles of the normalized vertical displacements at the bottom of the FG layer with a thickness h/a = 1.0, and the vertical profiles of the normalized vertical stresses along the z-axis of the FG layer with a thickness h/a = 2.0 from the current study were compared with the benchmark solutions [44] for different values of \bar{m}, and the results are shown in Fig.2(a) and 2(b), respectively. In the numerical implementation, the normalized thickness of each segmented FG layer was set to 0.1, where the underlying layer was assumed to be extremely large for simulating a homogeneous half-space. Fig.2 shows that both solutions yielded similar results in terms of both the displacement and stress profiles for all values of \bar{m}.

To validate the proposed solution scheme based on the effect of surface energy, a layered elastic medium under a uniform vertical surface load of radius a, as proposed by Tirapat et al. [8], was considered. In their work, a Si (100) upper layer with a finite thickness h/a = 1.0 overlaid the half-space of Al (111) without slip. Their bulk and surface elastic constants are listed in Tab.1 [5,40], and the normalized loading radius was \bar{a} = a/Λ = 10. The radial profiles of the normalized vertical surface displacements and vertical stresses along the interface are presented in Fig.3(a) and Fig.3(b), respectively, for different values of the residual surface stress at the interface {\tau }_2^s. To demonstrate the capability of the proposed solution scheme in managing the multilayered structure, the Si (100) layer was equally discretized into 10 layers of the same material properties, whereas the Al (111) half-space was discretized into 10 layers of sufficiently large thickness. To simulate the layered half-space presented by Tirapat et al. [8], the surface energy effects were only considered at the top surface and at the interface between Si (100) and Al (111), whereas they were disregarded at the interfaces between the discretized sublayers. Similarly, the comparison shown in Fig.3 indicates that the two solutions agreed well in terms of both the vertical surface displacements and vertical stresses.

3.2 Bi-material multilayered medium

In this subsection, the effects of surface energy and the size dependency characteristics of a multilayered nanoscale medium are presented. A Si/Al multilayered medium, which can be applied extensively in micro- and nanoelectric devices [45], was considered in the current study. The multilayered medium was composed of two distinct materials, i.e., Si (100) and Al (111) stacked alternately along the total thickness, H, and placed on a rigid substrate. It was subjected to a uniformly distributed vertical load applied over a circular area of radius a at its surface, as illustrated in Fig.4. The overall thickness of the multilayered medium was assumed to be 1 µm, and the thicknesses of the Si (100) and Al (111) layers were assumed to be equal, i.e., h_1 = h_2 = 0.2 nm. The bulk material properties of Si (100) and Al (111), and the top surface properties of Si (100) are listed in Tab.1. In addition, the surface material properties of Al (111) listed in Tab.1 were assumed for all the interfaces. The normalized loading radius, \overline{a} = a/\mathrm{\Lambda } =1.0, and the normalized thickness of each layer, h_1/\mathrm{\Lambda } = h_2/\mathrm{\Lambda } = 1.195, were used in the numerical calculation, where \mathrm{\Lambda } is the material length scale of the top layer, i.e., the Si (100) layer.

Fig.5(a) and Fig.5(b) show the radial profiles of the normalized vertical surface displacements and vertical stresses at depth \overline{\mathit{\text{z}}} = 0.1, respectively, of the bi-material layer under three different vertical loading cases, namely, a uniformly distributed vertical load, loading A, and loading B. All three loading cases were applied over a circular area of radius a with the same amount of total applied force. The functions of loading A and loading B were obtained from the contact pressure distribution under a rigid frictionless cylindrical and a paraboloidal indenter on the surface of an elastic half-space. They are expressed as follows, respectively [46].

The classical solutions for the vertical surface displacements due to loading A and loading B, as shown in Fig.5(a), show flat and paraboloid shapes, respectively. The contribution of the surface/interface energy effects is clearly presented in Fig.5(a); in particular, the bi-material layer became considerably stiffer owing to the surface stresses, which resulted in reduced vertical surface displacements in all loading cases. In addition, the numerical results presented in Fig.5(a) show that the maximum displacement occurred under the flat-ended loading case (loading A), and a kink occurred at the edge of the loading; meanwhile, the displacement was minimal under the paraboloid loading case (loading B). The surface energy effects were significant and resulted in reduced vertical stresses, as shown in Fig.5(b), particularly beneath the surface loading, i.e., r/a < 1.0.

To investigate the size-dependent characteristics of the multilayered medium, the normalized vertical surface displacements and vertical stresses at the depth of z = 0.1 nm of the Si/Al layered system under a uniformly distributed vertical surface loading were plotted against the size of the loading radius \overline{a}, as shown in Fig.6(a) and Fig.6(b), respectively. The solutions were plotted at r/a = 0.5, and the ratio H/a was maintained constant for all values of \overline{a}. Based on Fig.6, when surface energy was present, both the vertical displacements and stresses depended significantly on the size of the loading radius \overline{a}, whereas the classical solutions were size independent. The discrepancy between the present and classical solutions decreased gradually as \overline{a} increased and became insignificant when \overline{a} > 4.0; in other words, the effects of surface/interface energy are negligible for high value of \overline{a}.

3.3 Functionally graded layer

Next, the mechanical behavior of the FG layer was investigated by incorporating surface energy effects. Consider an FG layer of finite thickness h_1 completely attached to an underlying homogeneous elastic layer of thickness h_2, overlaying a rigid base, and subjected to a uniform vertical patch load of radius a at its top surface, as shown in Fig.7. The material properties of the FG layer were assumed to vary in the thickness direction, where the material properties of Si (100) and Al (111) were assumed based on depths of z = 0.0 and z = h_1 respectively. In this subsection, the overall normalized thickness H/a = 1.0 and normalized loading radius \overline{a} = 1.0 are used.

Fig.8(a) and Fig.8(b) show the radial profiles of the normalized vertical surface displacements and vertical stresses at the depth \overline{\mathit{\text{z}}} = 0.1 of an FG layer with different thickness ratios, i.e., h_2/h_1 = 1.5, 4.0, and 9.0. The elastic properties of the FG layer were assumed to reflect the exponential function L\left(\mathit{\text{z}}\right)=L_0{\mathrm{e}}^{\overline{m}\mathit{\text{z}}}, where \overline{m}={\displaystyle \frac{1}{h_1}}\mathrm{l}\mathrm{o}\mathrm{g}{\displaystyle \frac{L^{\ast }}{L_0}} is the grading constant, and L_0 and L^{\ast } represent the Lamé constants of Si (100) and Al (111), respectively. In the current study, the FG layer was segmented into 10 sublayers of equal thickness, with their elastic properties assigned a constant value obtained from the exponential function above at mid-height. In addition, the surface elastic properties of Si (100) and Al (111) were designated for the top surface and the interface between the FG layer and elastic substrate, whereas the effects of surface were disregarded at each interface of the FG sublayers.

The vertical surface displacements shown in Fig.8(a) indicate that both the classical and current solutions depend on the variation in the thickness ratio, particularly in the loading region, i.e., when r/a ≤ 1.0. Because the values of the elastic properties of Al (111) are lower than those of Si (100), the stiffness of this FG system is lower for a higher ratio of h_2/h_1 at a fixed value of H/a, thus resulting in a larger vertical surface displacement. By contrast, the variation in the thickness ratio does not significantly affect the vertical stresses of the FG layer, particularly in the classical case, as shown in Fig.8(b). In addition, the effects of surface energy reduce both the surface vertical displacements and vertical stresses under the loading region, and the discrepancy between the two solutions is negligible when r/a > 2.0.

For problems involving FG materials, the grading function is the key parameter that determines the mechanical behavior of those materials. To investigate the effect of the grading function, the radial profiles of the normalized vertical surface displacements and vertical stresses at the interface were obtained for different grading functions, namely, exponential, linear, and power-law distributions, as shown in Fig.9. The linearly distributed grading function is expressed as L\left(\mathit{\text{z}}\right)=L_0(1+\overline{m}\mathit{\text{z}}), where \overline{m}=(L^{\ast }-L_0)/\left(L_0{h}_1\right), and the power-law distribution is expressed in the form L\left(\mathit{\text{z}}\right)=L_0(1+\mathit{\text{z}}/h_1{)}^{\overline{m}}, where \overline{m}={\mathrm{l}\mathrm{o}\mathrm{g}}_2\left(L^{\ast }/L_0\right). Here, hypothetical materials with {\mu }_0 = 100 GPa, {\lambda }_0/{\mu }_0 = 1.5, {\mu }^{\ast } = 10 GPa, and {\lambda }^{\ast }/{\mu }^{\ast } = 1.5, along with the surface and interface material parameters of Si (100) and Al (111), respectively, were assumed. In addition, thickness values of \overline{h_1} = 0.4 and \overline{h_2} = 0.6 were considered in the numerical calculation.

Based on Fig.9, the effect of the grading function is more prominent on the vertical surface displacements than on the vertical stresses. The results of the vertical surface displacements shown in Fig.9(a) indicate that the FG layer with a linear distribution grading function is the stiffest layer, followed by the FG layers with exponential and power-law grading functions. However, results from extensive numerical experiments show that the effect of the grading function appeared only in extreme cases, where the difference between the bulk elastic properties at the top surface and interface was significant, as well as in cases involving extremely thick FG layers.

3.4 Frictionless indentation on bi-material multilayered media

Nanoindentations are used to determine the mechanical properties of extremely small systems such as nanocoatings, NEMS, and MEMS. In this subsection, the indentation problem is examined to exemplify the application of the proposed solution scheme. Fig.10 shows a rigid frictionless cylindrical indenter on a bi-material multilayered medium comprising Si (100) and Al (111). The shear surface traction on the contact surface vanishes, and the normal contact pressure can be computed by solving the boundary condition of the vertical surface displacement in the contact region, as expressed in the following equation.

(14){u_{\mathit{\text{z}}}|}_{\mathit{\text{z}}=0}=\underset{0}{\stackrel{a}{\int }}{U}_{\mathit{\text{z}}}^N\left(r,r^{\ast }\right)p\left(r^{\ast }\right)\mathrm{d}{r}^{\ast }=d,\text{0}\le r\le a,

where U_{\mathit{\text{z}}}^N\left(r,r^{\ast }\right) denotes the displacement Green’s function associated with the vertical surface displacement at point r under a unit normal ring load applied at radius r^{\ast } on the surface, p the contact pressure, and d the indentation depth.

To solve Eq. (14), the contact area is segmented into n annular ring elements, and the contact pressure corresponding to each element is assumed to be constant, as follows.

(15)p\left(r^{\ast }\right)=p_j,r_{j-1}^{\ast }\le r^{\ast }\le r_j^{\ast },j=1,2,...,n,

where p_j is the contact pressure at the middle of segmented element j.

By employing the collocation technique, Eq. (14) can be rewritten as

(16)\sum _{j=1}^n{U_{\mathit{\text{z}}}^{\mathrm{p}}\left(r_i,r_j^{\ast }\right)p}_j=d,i=1,2,...,n,

where U_{\mathit{\text{z}}}^{\mathrm{p}}\left(r_i,r_j^{\ast }\right) denotes the vertical surface displacement at r_i due to a uniformly distributed unit vertical surface load applied in the annular region r_{j-1}^{\ast }\le r^{\ast }\le r_j^{\ast }; it can be calculated by performing the inverse Hankel transform to the first component of the global displacement vector U_G (see Eq. (12c)). By solving Eq. (16), the distribution of the contact pressure can be obtained.

To validate the proposed technique, the normalized contact pressure and vertical surface displacements of a nanolayer with thickness h and λ/μ = 2.226, Λ = 1 nm, and τ^s = 5 N/m under a rigid cylindrical indenter of radius a obtained in the current study and by Intarit et al. [16] for various layer thicknesses were compared, as shown in Fig.11. In the present solutions, the layer was discretized into many sublayers of the same thickness (0.2 nm) and material properties. The normal ring load at r = a must be added to the formulation when considering surface energy effects to represent the ring induced by the displacement discontinuity at the punch edge. As shown in Fig.11, the current results are consistent with the benchmark solutions for all layer thicknesses.

Next, the effects of the surface/interface energy on the elastic responses of a multilayered system subjected to a rigid frictionless indenter was investigated. Consider a frictionless indentation of a Si/Al multilayer under a flat-ended cylindrical punch, as shown in Fig.10. All material properties of the layers, top surface, and interfaces were specified as presented in Subsection 3.2, except the residual surface stress of each interface {\tau }_2^{\mathrm{s}}, which was hypothetically assumed to be 0.5{\tau }_1^{\mathrm{s}}, 2.0{\tau }_1^{\mathrm{s}}, and 5.0{\tau }_1^{\mathrm{s}} to demonstrate the effects of interface energy. Here, {\tau }_1^{\mathrm{s}} is the residual surface stress of the top surface, i.e., the Si (100) layer. As shown in Fig.12(a), the normalized contact pressure increases with {\tau }_2^{\mathrm{s}}, and a singularity occurs in the contact pressure at the punch edge in both cases for all values of {\tau }_2^{\mathrm{s}}. Owing to the effects of interface energy in a multilayered system, which do not occur in a homogeneous elastic layer, the normalized contact pressure in this study was higher than that in the classical case [16]. In addition, Fig.12(b) shows that the normalized vertical surface displacement decreases as {\tau }_2^{\mathrm{s}} increases. This implies that the multilayered medium became stiffer, which necessitated a higher indentation force to deliver the same indentation depth for a higher value of {\tau }_2^{\mathrm{s}}.

4 Conclusions

A mathematical model for solving a multilayered nanoscale structure subjected to axisymmetric loading was proposed herein based on the G–M surface elasticity theory. An exact stiffness matrix approach was adopted along with Love’s strain function and the Hankel transform technique. Extensive parametric studies were conducted to demonstrate the capability of the solution scheme and the effects of the surface/interface energy on the elastic responses of the multilayered nanoscale medium. The numerical results showed that the surface/interface energy effects, particularly the contribution of the residual surface stress of both the top surface and interfaces, significantly affected the elastic fields and resulted in the size-dependent behavior of the multilayered nanoscale structures. The present solution can be employed as a benchmark for verifying numerical solutions based on the BEM and FEM to solve practical problems involving multilayered nanoscale structures.

References

[1] LacazeP C. Nanotechnologies: Concepts, Production and Applications. London: Wiley, 2012

[2] Wong E W, Sheehan P E, Lieber C M. Nanobeam mechanics: Elasticity, strength, and toughness of nanorods and nanotubes. Science, 1997, 277(5334): 1971–1975

[3] Gurtin M E, Ian Murdoch A. A Continuum theory of elastic material surfaces. Archive for Rational Mechanics and Analysis, 1975, 57(4): 291–323

[4] Gurtin M E, Ian Murdoch A. Surface stress in solids. International Journal of Solids and Structures, 1978, 14(6): 431–440

[5] Miller R E, Shenoy V B. Size-dependent elastic properties of nanosized structural elements. Nanotechnology, 2000, 11(3): 139–147

[6] Shenoy V B. Atomistic calculations of elastic properties of metallic FCC crystal surfaces. Physical Review B: Condensed Matter and Materials Physics, 2005, 71(9): 094104

[7] HeL HLimC WWuB S. A continuum model for size-dependent deformation of elastic films of nano-scale thickness. International Journal of Solids and Structures, 2004, 41(3−4): 847−857

[8] Tirapat S, Senjuntichai T, Rungamornrat J. Influence of surface energy effects on elastic fields of a layered elastic medium under surface loading. Advances in Materials Science and Engineering, 2017, 2017: 1–11

[9] Lu P, He L H, Lee H P, Lu C. Thin plate theory including surface effects. International Journal of Solids and Structures, 2006, 43(16): 4631–4647

[10] Sapsathiarn Y, Rajapakse R K N D. Finite-element modeling of circular nanoplates. Journal of Nanomechanics & Micromechanics, 2013, 3(3): 59–66

[11] Nanthakumar S S, Valizadeh N, Park H S, Rabczuk T. Surface effects on shape and topology optimization of nanostructures. Computational Mechanics, 2015, 56(1): 97–112

[12] Nguyen T B, Rungamornrat J, Senjuntichai T, Wijeyewickrema A C. FEM-SGBEM coupling for modeling of mode-I planar cracks in three-dimensional elastic media with residual surface tension effects. Engineering Analysis with Boundary Elements, 2015, 55: 40–51

[13] Intarit P, Senjuntichai T, Rungamornrat J, Rajapakse R K N D. Penny-shaped crack in elastic medium with surface energy effects. Acta Mechanica, 2017, 228(2): 617–630

[14] Tian L, Rajapakse R K N D. Analytical Solution for size-dependent elastic field of a nanoscale circular inhomogeneity. Journal of Applied Mechanics, 2007, 74(3): 568–574

[15] Pinyochotiwong Y, Rungamornrat J, Senjuntichai T. Rigid frictionless indentation on elastic half space with influence of surface stresses. International Journal of Engineering Science, 2013, 71: 15–35

[16] Intarit P, Senjuntichai T, Rungamornrat J. Elastic layer under axisymmetric indentation and surface energy effects. Zeitschrift für Angewandte Mathematik und Physik, 2018, 69: 29

[17] Tirapat S, Senjuntichai T, Rungamornrat J, Rajapakse R K N D. Indentation of a nanolayer on a substrate by a rigid cylinder in adhesive contact. Acta Mechanica, 2020, 231(8): 3235–3246

[18] Intarit P, Senjuntichai T, Rungamornrat J, Limkatanyu S. Influence of frictional contact on indentation of elastic layer under surface energy effects. Mechanics Research Communications, 2020, 110: 103622

[19] Eringen A C. Nonlocal polar elastic continua. International Journal of Engineering Science, 1972, 10(1): 1–16

[20] Eringen A C, Edelen D G B. On nonlocal elasticity. International Journal of Engineering Science, 1972, 10(3): 233–248

[21] Mindlin R D, Tiersten H F. Effects of couple-stresses in linear elasticity. Archive for Rational Mechanics and Analysis, 1962, 11(1): 415–448

[22] Toupin R A. Theory of elasticity with couple-stress. Archive for Rational Mechanics and Analysis, 1964, 17(2): 85–112

[23] Mindlin R D. Micro-structure in linear elasticity. Archive for Rational Mechanics and Analysis, 1964, 16(1): 51–78

[24] Wongviboonsin W, Gourgiotis P A, Van C N, Limkatanyu S, Rungamornrat J. Size effects in two-dimensional layered materials modeled by couple stress elasticity. Frontiers of Structural and Civil Engineering, 2021, 15(2): 425–443

[25] Samaniego E, Anitescu C, Goswami S, Nguyen-Thanh V M, Guo H, Hamdia K, Zhuang X, Rabczuk T. An Energy Approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications. Computer Methods in Applied Mechanics and Engineering, 2020, 362: 112790

[26] Anitescu C, Atroshchenko E, Alajlan N, Rabczuk T. Artificial neural network methods for the solution of second order boundary value problems. Computers, Materials & Continua, 2019, 59(1): 345–359

[27] Nguyen-Thanh V M, Anitescu C, Alajlan N, Rabczuk T, Zhuang X. Parametric deep energy approach for elasticity accounting for strain gradient effects. Computer Methods in Applied Mechanics and Engineering, 2021, 386: 114096

[28] GürdalZHaftkaR THajelaP. Design and Optimization of Laminated Composite Materials. New York: Wiley, 1999

[29] Senjuntichai T, Sornpakdee N, Teerawong J, Rajapakse R K N D. Time dependent response of an axially loaded elastic bar in a multi-layered poroelastic medium. Journal of Engineering Mechanics, 2007, 133(5): 578–587

[30] Senjuntichai T, Kaewjuea W. Dynamic response of multiple flexible strips on a multilayered poroelastic half-plane. Journal of Mechanics of Materials and Structures, 2008, 3(10): 1885–1901

[31] Ba Z N, Kang Z Q, Lee V W. Plane strain dynamic responses of a multi-layered transversely isotropic saturated half-space. International Journal of Engineering Science, 2017, 119: 55–77

[32] Senjuntichai T, Keawsawasvong S, Plangmal R. Vertical vibrations of rigid foundations of arbitrary shape in a multi-layered poroelastic medium. Computers and Geotechnics, 2018, 100: 121–134

[33] Ba Z N, Sang Q Z, Liang J W. Seismic analysis of a lined tunnel in a multi-layered TI saturated half-space due to qP1-and qSV-waves. Tunnelling and Underground Space Technology, 2022, 119: 104248

[34] ZhuX YLiuX JZongR LZengFPanF. Microstructure and mechanical properties of nanoscale Cu/Ni multilayers. Materials Science and Engineering A, 2010, 527(4−5): 1243−1248

[35] Kern W, Schuegraf K K. Handbook of Thin Film Deposition Processes and Techniques. 2nd ed. Norwich: William Andrew Publishing, 2001, 11–43

[36] Senjuntichai T, Rajapakse R K N D. Exact stiffness method for quasi-statics of a multi-layered poroelastic medium. International Journal of Solids and Structures, 1995, 32(11): 1535–1553

[37] SneddonI N. Fourier Transform. New York: McGraw-Hill, 1951

[38] Povstenko Y Z. Theoretical investigation of phenomena caused by heterogeneous surface tension in solids. Journal of the Mechanics and Physics of Solids, 1993, 41(9): 1499–1514

[39] Wang Z Q, Zhao Y P, Huang Z P. The effects of surface tension on the elastic properties of nano structures. International Journal of Engineering Science, 2010, 48(2): 140–150

[40] MeyersM AChawlaK K. Mechanical Behavior of Materials. 2nd ed. Cambridge: Cambridge University Press, 2008

[41] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

[42] Hamdia K, Ghasemi H, Zhuang X, Alajlan N, Rabczuk T. Sensitivity and uncertainty analysis for flexoelectric nanostructures. Computer Methods in Applied Mechanics and Engineering, 2018, 337: 95–109

[43] PiessensRDoncker-KapengaEUberhuberC WKahanerD K. QUADPACK: A Subroutine Package for Automatic Integration. Berlin: Springer, 1983

[44] Katebi A, Selvadurai A P S. Undrained behaviour of a non-homogeneous elastic medium: The influence of variations in the elastic shear modulus with depth. Geotechnique, 2013, 63(13): 1159–1169

[45] NakayamaTYamamotoKSatohHKonnoT JClemensB MSinclairR. Structure and corrosion properties of Al/Si and Fe/Zr multilayers. Materials Science and Engineering A, 1995, 198(1−2): 19−24

[46] Sneddon I N. The relation between load and penetration in the axisymmetric Boussinesq problem for a punch of arbitrary profile. International Journal of Engineering Science, 1965, 3(1): 47–57

RIGHTS & PERMISSIONS

Higher Education Press