Abstract
Crack initiation and propagation in a nanostructured nickel film were studied by molecular dynamic simulation as well as an interatomic-potential-based continuum approach. In the molecular dynamic simulation, the interatomic potential was described by using Embedded Atom Method (EAM), and a reduced 2D plane model was employed to simulate the mechanical behavior of nanofilms. Atomistic simulation shows that the reduced plane model in this paper can not only reveal the physical nature of crack initiation clearly but also give the critical time of crack initiation accurately as the continuum fracture analysis does. The normal stress and average atom energy at the crack tip which resulted from atomistic simulation at the time of crack initiation agree well with the analytical results. On the other hand, the crack propagation in nanofilms was studied by interatomic-potential-based continuum fracture mechanics analysis based on Griffith criterion. The coupled continuum-atomic analysis can predict the crack initiation and atomic stress accurately. Continuum analysis with material property parameters determined by interatomic potential is proved to be another promising way of solving failure problem on nanoscale.
1. Introduction
Nanofilms have been more and more widely used in various industries up to date. An open and active challenge in the mechanical behavior of nanofilms is the mechanism of crack initiation and propagation. Recently, a wide variety of methods have been employed to measure the crack initiation, propagation, or healing in nanofilms [1–3] and analyze the mechanism of crack initiation [4–6] and the fracture behavior of nanofilms [7–9]. On the other hand, continuum-based theories of fracture mechanics have provided plenty of energy [10] and force [11] criteria for computing the initiation and propagation process of observable fracture. Despite their valuable contributions to the development of fracture mechanics, these theories and criteria on macroscopic scale cannot account for more and more phenomena experimentally observed, especially crack characteristics on microscopic scale in the latest two decades. On the other hand, molecular dynamics (MD) method proposes an ab initio simulation for mechanical behavior of nanoscaled materials and structures, provided that the interatomic potential can be accurately described. With the development of computer hardware and emergence of smart algorithms, atomistic simulations of mechanical behavior and failure process of nanoscaled materials are performed with larger and larger amount of atoms. However, as the coarse view of physical phenomena and constitutive assumptions in continuum fracture mechanics, atomistic methods are still limited by length and time scales, which leads to a gap between continuum and atomistic analysis and leaves only a few problems adapted to both methods and the results can be compared directly. Crack propagation in nanometer size crystals is one of such problems [12] where the system is large enough for continuum analysis and the length and time scales are suitable for molecular dynamics simulations. From this point of view, making a multiscaled model of fracture or performing both continuum and atomistic analyses simultaneously has become a promising method to reveal the mechanical and physical nature of failure process of materials. Yang et al. [13] and Tan and Yang [14, 15] have put forward an idea of critical stress intensity factor rate to explain the tough-brittle transition during fracture process zone at nanoscale and made a combined atomistic/dislocation/continuum model to simulate the fracture of materials with emphasis on establishing compatible boundaries across different scales. In addition, Ozmat et al. [19] observed crack bifurcation and appearance of daughter cracks in creep experiments.
Gao et al. [12] have attempted to establish a linkage between continuum and atomistic analyses by applying dynamic elasticity methods to solve crack propagation problems with material properties determined by interatomic potentials used in atomistic simulation. By using this method the emergence of daughter crack at a distance ahead of the mother crack in atomistic simulation is captured, and the predictions based on continuum analysis agree remarkably well with atomistic results. Ma and Yang [17] investigated the fast deformation of nanocrystals and order-disorder transition. What is more, Rafii-Tabar et al. [18, 19] proposed an atomistic fracture model with two-dimensional triangular lattices to represent three-dimensional face-centered cubic (fcc) metals so as to reduce the large number of degrees of freedom in molecular dynamics simulation of crack propagation. In Rafii-Tabar’s work, the convergence, accuracy, and computational efficiency of employing a simplified 2D model to simulate 3D phenomena have been validated.
In the present study, a joint continuum and atomistic investigation of intersonic crack propagation in nanofilms is proposed, and a study of propagation of mode crack in a two-dimensional nanocrystalline nickel film is conducted using both molecular dynamics method and interatomic-potential-based continuum analysis. The problem selected is an edge-notched plane-stress nanocrystalline nickel film consisting of reduced two-dimensional triangular lattices with interatomic potential described by Embedded Atom Method [20]. In molecular dynamic simulation, a fixed tensile strain rate was imposed on atoms along the direction normal to the preexisting crack, until the critical time when the crack becomes instable and starts to propagate. Sometime after the critical time, if the external strain rate keeps unchanged, the bifurcation and quadratic bifurcation of crack in front of the original crack tip will appear, just as Ozmat et al. [16] have reported. On the other hand, the dynamic elasticity methods introduced by Gao et al. [12] were employed to solve the plane-stress field near the crack tip and obtain the critical time for crack initiation with the same geometry and loading conditions as the atomistic simulation. The material property parameters used in continuum analysis resulted from the interatomic potential. Such a linkage between continuum and atomistic analysis can validate the agreement of atomistic simulation results with traditional dynamic elasticity analysis.
2. Molecular Dynamics Approach
Figure 1 shows the molecular dynamics simulation model. It is a two-dimensional slab with 60 atoms along the horizontal length defining the direction and 100 atoms along the vertical length as the direction. A horizontal notch of 11-atom length is cut midway from the left vertical boundary. Atoms right above the notch do not interact with the atoms under the notch in the model. The slab is composed of two-dimensional triangular lattices, which describe the atom arrangement of three-dimensional fcc metal, nickel, for the sake of simplification of computation. The model is set to be periodic along vertical orientation by applying periodic boundary condition (PBC), and the left and right boundaries of the model are set to be free surfaces. Interatomic potential of nickel is described by EAM [20, 21], with geometrical revision in accordance with the two-dimensional lattice. System atoms are assumed to interact with each other through the potential function where is the pair potential between atoms and and is the embedding energy of atom into electron density , and they are presented as follows: Here, stand for the contribution of atom at the distance from atom to the electron density . Constants , , , , and are determined by material physical prosperities including the Born stability, cohesive energy, elastic constants, formation energy of a vacancy, and stacking fault energy. Constants and are smoothly cut-off distance of and , respectively, and in this paper are both set to be , which is between the third nearest neighbor atoms and the fourth nearest neighbor atoms in the two-dimensional model. Constant is the two-dimensional triangular lattice parameter of nickel.
The motion of each atom is predicted through Newton’s law: and the numerical integrating method chosen here is the velocity form of leapfrog Verlet arithmetic [22]: which enhances the computation efficiency with only variables at one time stored in the memory once.
Values throughout this paper are expressed in terms of reduced units as follows: time is scaled by the selective time step, 3 fs, which is proved to be the optimum time step size through stability and efficiency tests, energies are scaled by the minimum value of EAM potential, 4.227 × 10−12 erg, mass is scaled by the atomic mass, and lengths are scaled by two-dimensional triangular lattice parameter , 2.489 angstroms ( times of three-dimensional nickel fcc lattice parameter).
The initial temperature of the model was 0 K, which kept invariant during the simulation to avoid atomic heat activation and did not allow for the appearance of plastic zone and brittle-to-ductile transition during fracture simulation. A static tensile strain rate of 4 × 109/s along direction was imposed on the atoms of the outermost rows; an elastic atomic model of semi-infinite slab with a boundary mode crack then came into being finally. The corresponding atomic stress and energy (kinetic energy, potential energy, and total energy), atom arrangement along the crack surface and near the fracture field, and the form change of crack tip are recorded instantaneously and written into output files every 100 time steps.
3. Continuum Analysis and Atomistic Simulation Results
Gao [23] deducted the moduli constants of two-dimensional triangular lattice. According to the method proposed in Gao’s work, the elastic constants from the EAM interatomic potential used in atomistic simulation can be obtained as follows: where , , and are Young’s modulus, shear modulus, and Poisson ratio, respectively. The density of triangular lattice is then, we can get the longitudinal and shear wave speed as The fracture surface energy can be expressed as the stored energy in the breaking bonds as crack grows. In this paper, the cut-off distance of interatomic potential was 2.25, as mentioned above. So we snap two nearest atomic bonds and two second nearest atomic bonds to calculate the fracture surface energy per atom: According to elastic fracture mechanics theories, the stress intensity factor and energy release rate for a stationary plane-stress mode crack are With Griffith criterion on the relationship between crack tip energy release rate and crack surface energy at the time that crack starts to propagate, , we obtain the critical time for mode crack initiation: From (6), (8), (9), and (12), the given tensile strain rate above, and linear elastic constitutive relation, the theoretical critical time for crack initiation can be calculated as On the other hand, we can observe the variation of atom arrangement in fracture field during the loading period, compare the atomistic simulation results with continuum analysis results, and test the accuracy of our two-dimensional fracture model on the prediction of crack initiation.
Figure 2 shows the atom arrangement every 2000 steps during the atomistic simulation. At step 2000, atoms in the model are still in regular arrangement as in the initial time. When , we find that atoms along the crack surfaces start to depart from the initial lattice location, and there is a triangular field around the crack, in which the atoms are not in their perfect lattice location anymore. The triangular field is symmetrical with the crack extending line as the centerline. When , the triangular field enlarges (as shown in Figure 2(c)) and, at step 8000, two new vacant atom sites appear symmetrically at the crack tip, which indicates that the crack starts to propagate. It is shown that the atomistic simulation results agree well with the continuum dynamic elasticity analysis on the prediction of critical time for mode crack initiation. What is more, atomistic simulation reveals the physical nature of crack propagation compared with continuum analysis.
(a)
(b)
(c)
(d)
With the stress intensity factor , we can obtain the stress field near the crack tip through linear elastic fracture mechanics analysis. Subjected to constant remote normal traction rate, the tensile stress normal to the crack direction near the crack tip is then, with (10), the normal tensile stress at the crack tip is given by In addition, with Griffith fracture criterion on brittle crack initiation, the average fracture surface energy per atom when the crack starts to grow is Finally, the substitution of given tensile stress rate and calculated constants above into (15) and (16) brings the functions to determine the variation of normal tensile stress and fracture surface energy per atom at the crack tip as Thus, we can work out the tensile stress normal to the crack and the average fracture surface energy in the crack tip at the time of crack initiation (when according to the theoretical analysis above). Assuming that , that is, the location at which we calculate the stress is one-atom distance (the nearest neighbor distance in the two-dimensional lattice) from the calculating coordinate center of the crack tip (in fact it is just one-atom distance from which the crack initiation occurs, as shown in Figure 2(d)), then the normal tensile stress is obtained as and, at the same time, the energy stored per atom on average in the breaking bonds is
Figure 3 shows the atomistic simulation results of stress and energy. Figure 3(a) describes the normal stress variation of the crack tip atoms, and Figure 3(b) tells the atom total energy variation of the same atoms. According to the atomistic simulation output files, when , the average normal tensile stress and atom total energy of the crack tip atoms are and , respectively. Compared with the data in (18) and (19), the atomistic simulation results appear to agree well with continuum analysis results again.
(a) Variation of crack tip atomic stress
(b) Variation of atom total energy
From continuum analysis and atomistic simulation and the comparison of their results, we can conclude that the two-dimensional molecular dynamics model of nanonickel can predict the time of crack initiation accurately and close to the analysis result based on continuum fracture mechanics theories. And the results of tensile stress normal to the crack and the atom total energy of the crack tip atoms at the time of crack initiation by two different methods are also in good agreement with each other.
After a period of crack propagation, a symmetrical bifurcation ahead of the crack tip appears at , as shown in Figure 4. From the instantly caught pictures it can be found that the crack propagates not only smoothly along the crack tip but also along normal direction. The position where the daughter crack was born is 12-atom distance ahead of the mother crack tip in direction and 5-atom distance from the horizontally extended line of the mother crack in direction. The symmetric daughter cracks are along the direction of departure from the horizontal direction. The daughter cracks in quadratic bifurcation are almost parallel to the cracks in the first bifurcation and two incontinous cracks in each side of the bifurcation. There are also other atomic cavities in the region near the bifurcation crack tip, which accords with previously reported experimental observations [16].
4. Conclusions
With molecular dynamics method and continuum analysis, we have studied the mode crack initiation in two-dimensional nanonickel slab. The atomistic simulation reveals the sequent events of atoms around the crack and the simulation results agree well with the prediction of continuum fracture mechanics analysis based on Griffith criterion. In atomistic simulation, the critical time for crack initiation, the atomic normal tensile stress, and atom total energy at the crack tip when crack starts to propagate are , , and , while the data are , , and , respectively, from the continuum analysis. The remarkable agreement of results from two different methods indicates the validity of the mode crack simulation model which describes the interatomic actions by Embedded Atom Method and nickel crystal configuration by reduced two-dimesnional triangular lattices, with parameters modified accordingly. On the other hand, continuum fracture mechanics analysis with physical property parameters determined by interatomic potenial has been proved to be another effective way to solve nanoscale material fracture problems. Bifucation of crack ahead of crack tip appears, which agrees well with experimental observations. The coupling method of continuum analysis and atomistic simulation is expected to be a promising way of studying material failure process, from physical and theoretical point of view simultaneously.
Conflict of Interests
The authors do not have any direct financial relation with the commercial identity that might lead to a conflict of interests for any of them.
Acknowledgment
This work was partially supported by the National Natural Science Foundation of China through Grants no. 11132003 and no. 51209079.