Abstract
This paper reports a systematic study on the elastic property of bulk silicon nanomaterials using the atomic finite element method. The Tersoff-Brenner potential is used to describe the interaction between silicon atoms, and the atomic finite element method is constructed in a computational scheme similar to the continuum finite element method. Young’s modulus and Poisson ratio are calculated for 100], 110], and 111] silicon nanowires that are treated as three-dimensional structures. It is found that the nanowire possesses the lowest Young’s modulus along the 100] direction, while the 110] nanowire has the highest value with the same radius. The bending deformation of 100] silicon nanowire is also modeled, and the bending stiffness is calculated.
1. Introduction
Over the past decade, one-dimensional nanostructures such as nanowires have been extensively studied through both theoretical and experimental methods. The most important one is the silicon nanowire (SiNW) that has been successfully applied in the nanoelectromechanical devices such as field effect transistors (FETs) [1–3] due to their unique optoelectronic and mechanical properties [4, 5]. Researchers have made significant contributions in the study of the elastic property of SiNWs, in which they are considered with different diameters and different axial or surface orientations. It is generally accepted that Young’s modulus tends to increase with the increasing diameter [6, 7].
The molecular dynamics and density functional theory [6–9] are the common atomic scale simulating methods to study the elastic property of nanostructures. However, their computational cost is very huge and they are valid only for the small size structures. In order to overcome this shortcoming, Liu et al. [10, 11] proposed an order-N atomic finite element method (AFEM) which can achieve the same accuracy with molecular mechanics but is much faster than the latter one. Sun and Tao [12] adopted this method to investigative the elastic properties of carbon, boron nitride, and silicon carbide nanotubes, and Young’s moduli are accurately obtained and the buckling behavior is better displayed. The present work extends its application to three-dimensional nanostructures in order to study the fine elastic property of SiNWs. An appropriate type of Tersoff-Brenner potential is employed to describe the atomic interaction. An atomic element contains 17 atoms for the bulk silicon crystal, and the global stiffness matrix and nonequilibrium force vector are assembled similar to the continuum finite element method (FEM). The nonlinear iteration is implemented to determine the equilibrium.
2. Atomic Scale Modeling Method for Bulk Silicon
Bulk silicon has the diamond lattice structure as illustrated in Figure 1, and its elastic property has a great dependence on their growth directions. The present work focuses mainly on SiNWs grown along the 100], 110], and 111] crystallographic orientations, and Figure 2 shows the representative cross section and side view of the SiNWs with different growth directions.
![(a) 100]](https://static.hindawi.com/articles/jnm/volume-2016/8231592/figures/8231592.fig.002a.jpg)
(a) 100]
![(b) 110]](https://static.hindawi.com/articles/jnm/volume-2016/8231592/figures/8231592.fig.002b.jpg)
(b) 110]
![(c) 111]](https://static.hindawi.com/articles/jnm/volume-2016/8231592/figures/8231592.fig.002c.jpg)
(c) 111]
In order to describe the Si-Si interaction, the Tersoff-Brenner potential is adopted with the parameter set refined by Erhart and Albe [13]. Tersoff-Brenner potential is a multibody potential which can be used to describe the interaction between atoms including C, B, H, Si, and N: where and are the repulsive pair potential and attractive pair potential, respectively; is the distance from atom to atom ; is the bond order function. All above functions have analytic forms where is the cutoff functionThe angular function is given byThe potential parameters are listed in Table 1.
The basic idea of AFEM is to treat atoms as the continuum FEM nodes, and each element is characterized by a set of discrete atoms [10, 11]. For a system of atoms, the energy stored in the atomic bonds can be evaluated using Tersoff-Brenner potential, and the atomic positions are determined by minimizing the system energy in which , is the position of atom , and is the external force exerted on atom .
Giving Taylor expansion of and substituting it intowe havewhere is displacement increment, and are, respectively, the stiffness matrix and nonequilibrium force vector given by
For the present nonlinear system, (7) is solved iteratively until reaches zero. Here, Newton method is used to solve the problem.
In AFEM, the choice of element depends on the atomic structure and nature of atomistic interaction. In the studies for carbon nanotubes by Liu et al. [10, 11], an AFEM element containing 10 atoms is developed in order to capture the multibody interaction. For the current diamond lattice structure and multibody interaction, an AFEM element should contain 17 atoms that include 4 nearest-neighbor atoms (red balls) and 12 second nearest-neighbor atoms (blue balls) as shown in Figure 3. Therefore, the element stiffness matrix and the nonequilibrium force vector are given by
In the present study, all work, including assembling the stiffness matrix and force vector and solving the equation system, is performed with our Fortran codes. The following steps are used. Firstly, the initial configuration of SiNWs with a uniform bond length is constructed, and the obtained coordinates and bond information are stored in an array. Using the obtained array, the information about the first and second neighbor atoms for each atom can be found and is stored in another array. Secondly, boundary conditions are applied to SiNWs to get the initial equilibrium coordinates of the system. Here, we restrain one side of SiNWs and make the other side free until the system returns to the equilibrium configuration. Table 2 gives the bond lengths of different types of SiNWs when the potential parameters in [13] are employed. These results are similar to the work in [14]. In the third step, displacement field is applied to the initial equilibrium coordinates in order to simulate the axial tension or compression. It is achieved by completely fixing one side of SiNWs and incrementally imposing an axial movement at the other end. The length of the nanowire is changed by 0.01 nm per loading step. In the final step, the total energy of SiNWs is obtained, and it is treated as a function of strain. With the first- and second-order derivatives, Newton iteration method is thus applied to obtain the equilibrium state, and 3–5 iterative steps can generally achieve a good convergence. Therefore, the present method is far faster than molecular dynamic method and is very efficient for the nanostructures with a larger number of atoms.
3. Results and Discussions
In continuum mechanics, the material can be represented by two independent constants, namely, Young’s modulus and Poisson’s ratio . For a material undergoing a uniaxial deformation, is defined as [11, 15]where is the volume corresponding to the initial equilibrium configurations of SiNWs; is the total strain energy under tensile or compressive deformation, and is strain.
The Poisson ratio is calculated as where and are the mean diameters of SiNWs corresponding to the initial and deformed equilibrium configurations, respectively. They can be derived from , where is the cross-sectional area.
Young’s modulus for the bulk silicon crystal under tension along 100] directions is first calculated in this paper, and it has the bulk structure with the side length 2.16 nm. The total potential energy is plotted versus strain in Figure 4. Using the polynomial curve fitting, the obtained is 134.03 GPa. For 110] and 111] directions, the present computational results are, respectively, 190.2 GPa and 178.3 GPa, which drop in the literature reports. Ma et al. [6] used all-electron density functional theory to calculate the binding energy, heat of formation, and Young’s modulus of hydrogen-passivated SiNWs with various diameters and orientations, and Young’s moduli of bulk silicon for the 100], 110], and 111] directions are 133.24 GPa, 157.66 GPa, and 163.58 GPa, respectively. Lee and Rudd [15] found that the bulk value of Young’s modulus is 122.53 GPa for the 100] direction by the density functional theory.
Next, the numerical computation is carried out for Young’s modulus and Poisson ratio of SiNWs with different dimensions. Table 3 shows the Poisson ratio for 100], 110], and 111] SiNWs, which have a small dependence on diameter. The calculated Young’s moduli are plotted in Figure 5. Leu et al. [5] have given a systematic study of the mechanical property of strained small diameter silicon nanowires using ab initio density functional theory, and the values of Young’s modulus are 139–153 GPa and 43–131 GPa for 110] and 111] SiNWs, respectively. In the researches by Ma et al. [6] and Lee and Rudd [15], Young’s modulus increases while the wire mean diameter increases, and at certain critical size it approaches the bulk value. The present results agree well with these previous studies. It can also be seen from Figure 5 that, with the same mean diameter, nanowires with different grown directions have distinctly different Young’s moduli, and the wires along 110] direction have the highest value while 100] wires have the lowest value.
Besides the axial strain, the simulation has also been carried out for the bending of SiNWs. The bending stiffness of a SiNW can be calculated aswhere is the length and is the curvature of the bent SiNWs and is related to the bending angle as [16].
The bending deformation can be achieved by incrementally rotating the two end planes in the opposite directions with respect to a line that is perpendicular to the axis of the undeformed SiNW and passing through its center. During the simulating process, the two end planes are reversely rotated 0.1 rad per loading step. Two 100] SiNWs with different diameters are studied in this paper. The change of the total energy is plotted as a function of the bending angle in Figure 6. Using polynomial curve fitting, we can obtainThe bending stiffness can be calculated by replacing with . The bending stiffness is, respectively, 5432.7 and 1964.2 ev·nm for SiNWs with the diameters 2.43 nm and 1.83 nm. Menon et al. [17] investigated the bending stiffness of T and C type nanowires, and their computed results are 12.2 × 10−6 ev·m and 8.5 × 10−6 ev·m, respectively.
4. Conclusions
The AFEM element is developed for the silicon nanowire, and a basic element contains 17 atoms while it is treated as the bulk structure. Young’s modulus and Poisson ratio are calculated for 100], 110], and 111] SiNWs, and the bending stiffness is also obtained. AFEM is proved to be an accurate atomic simulation method for the bulk nanomaterials, and it is found that the 110] nanowire has the largest Young’s modulus while the 100] nanowire has the smallest value. Newton iteration method is applied to obtain the equilibrium state, in which the first- and second-order derivatives are used and 3–5 iterative steps can achieve a good convergence, and the developed method is much faster than the molecular dynamic. In order to take advantage of the order-N nature of AFEM, further studies can be made to combine AFEM and the classical FEM software to solve problem with a large number of degrees, for example, the nanoindentation of bulk monocrystalline silicon.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
The work in this research was supported by the Natural Science Foundation of China (Grant no. 11472316), Innovation Scientists and Technicians Troop Construction Projects of Henan Province (Project no. 164200510020), and Natural Science Foundation from Department of Education of Henan Province (Project no. 13B560304).