Metallic nanolayered composites (MNCs) are nanostructured materials possessing impressive mechanical properties such as high strength, hardness and fatigue resistance, which can be fabricated by deposition or severe plastic deformation processes [1][2][3][4][5][6] . MNCs are composed of alternating layers of two or more metallic phases and the layer thickness of each phase are generally less than 100 nm. The interfaces between each layer play an critical role on the deformation of MNCs, as they can act as barriers, sink and sources of dislocations and vacancies [7][8][9] . Both experimental and computational studies have illustrated the remarkable thicknessdependent strength in MNCs that can be predicted by the confined layer slip (CLS) model [10][11][12][13] . Similar to nanocrystalline (NC) metallic materials, MNCs also exhibits an inverse relationship between the strength and elongation 12,14 . To explore the fracture mechanisms in MNCs, Zhu at el. 15 examined the deformation zone ahead of the crack tip in the Cu/Ta MNCs and revealed a critical layer thickness, below which the fracture mode of the MNCs tends to be shearing failure.
Zhang at el. 16 studied the fracture behavior of Cu/Nb and Cu/Zr. Their experiment results demonstrated that as the layer thickness of Cu layer decreased below 60 nm, the fracture mode in MNCs transited from brittle opening fracture to shear fracture. Based on their experiment results, they claimed that the transition of fracture modes is dominated by the constraint of the soft Cu layer on the brittle Nb or Zr layer. Hattar et al. 17 demonstrated four fracture steps (crack deviation, layer necking, microvoid formation and crack blunting) during the crack propagation in Cu/Nb by using the in-situ transmission electron microscopy testing. Radchenko et al. 18 found that the crack propagation under three-point bending in 63 nm Cu/Nb MNCs was inhibited by the interface shear. Atomistic simulations performed by Li et al. 19 also illustrated the weak interface in Cu/Nb MNCs can hinder the crack propagation in Nb phase by interface shear and dislocation nucleation from interface.
Based on previous studies, we can see that in most MNCs, the plastic deformation ability is limited by the thickness of the each single layered phase. The interface in the MNCs can not only influence the strength of MNCs but also affect their ductility and the fracture behaviors. However, it is still unclear how the interface in MNCs affects the fracture mode of MNCs and how the crack interacts with the interface, especially for the cracks cutting through the layers.
Atomistic simulations, such as molecular dynamics (MD) simulations, can shed light onto the failure mechanisms of NC materials by directly revealing the underlying atomic scale processes of the deformation and fracture. During the fracture, the crack propagation speed is comparable to speed of sound, and is well captured in the time frame of MD simulations. MD simulations have already been successfully used to explore the fracture behavior for NC metals, nanotwinned metals and gradient metals [20][21][22][23][24] . In this study, we perform MD simulations of the deformation in Cu/Nb MNCs with a preexisting crack to explore the crack propagation and fracture behavior of MNCs. Our simulation results show that the interface can increase the fracture resistance of the brittle Nb layer by providing the dislocation sources to generate the plastic strain at the front of the crack. The fracture resistance of both layers would be enhanced by increasing the layer thickness. Introducing grain boundaries (GBs) would decrease the fracture resistance of Nb layer, as the GBs act as the weakest location promoting the intergranular fracture.
To study the effect of interface and the coupled effect of interface and grain boundaries on fracture behavior of MNCs, we adopt multiple types of samples in our model which include two-phase samples of single crystalline (SX) Cu/Nb, polycrystalline (PX) Cu/Nb, and single phase samples of SX Cu, SX Nb, PX Cu, PX Nb. Figure 1 presents the examples of samples used in our calculations. Periodic boundary conditions have been applied along the thickness direction, which normally contains a large number of thin layers. As our study focuses on the mechanism of crack propagation within the nanolayers, this model may restrict a possibility of the shear localization that may take place at the final/catastrophic fracture stage. Fixed boundary conditions were set in other two directions under external loading. A pre-existing crack is created with the crack tip end at the center of each sample. To mimic the microstructure of Cu/Nb composites synthesized by the physical vapor deposition method (PVD) 25 , the crystallographic orientation for Cu layer is set as x- [11-2], y-[1-10], and z-[111], while the crystallographic orientation for Nb layer is set as x- [11-2], y-[111], and z- [1-10]. This configuration of Cu and Nb phases follows the Kurdjumov-Sachs (KS) orientation relationship 26 . The dimension of the simulation box was optimized to minimize the internal stresses. In the PX Cu/Nb samples, a modified Voronoi method 27 was used to create the nano grains with a hexagonal shape shown in Figure 1(b). The grain size for PX samples in this study set to be 10nm. We fixed one pair of Cu/Nb grains with the initial crystallographic orientation described above and rotated the neighboring grain pairs by 30°, 60° and 90° degree along the thickness direction shown in Figure 1 (b). In this way, stable large angle GBs can be created and the KS orientation relationship can still be maintained in each pairs of Cu/Nb grains, this method has also been used in previous MD simulations of polycrystalline Cu/Nb and Cu/Ag MNCs 12,28 . Finally, to explore how the layer thickness affects the fracture behavior of MNCs, we vary the layer thickness of both SX Cu/Nb and PX Cu/Nb samples from 5 nm to 20 nm in our model. the final propagation distance for Nb layer in 20nm SX Cu/Nb decrease from 21 to 6 nm. Thus, the thicker the sample, the better the crack resistance is in the Nb layer. This trend is also consistent with the results from previous experiment study by Kavarana et al. 14 , which demonstrated that that the ductility of MNCs increases with the bilayer thickness. The evolution of the Green strain in each case is consistent well with the crack propagation distance shown in Figure 3(a) and SIFs in Figure 3(b). The phase with higher SIF and shorter crack propagation distance normally contains larger Green strain. In metallic materials, there are two ways to release the excess elastic energy stored within the materials: i) creating new surfaces by opening cracks, ii) changing the shape of the material via plastic deformation. In the SX samples, the plastic deformation was mainly induced by the dislocation activities. Thus, larger Green strain under the same applied strain indicates more plastic deformation induced by dislocation activates that suppressed the crack propagation to create new surfaces. In the Cu/Nb samples, the interface provides dislocation sources for nucleating interfacial dislocations into each phase as shown in Figure 4 (c) and (d). Thus, the Nb layers in Cu/Nb samples own more Green strain than the SX Nb sample. In addition, thicker Cu/Nb samples carry larger Green stains in both Cu and Nb layers. That is because the confined layer slip (CLS) is the main dislocation activity in the SX Cu/Nb samples as shown in Figure 4 (d). The critical stress for activating confined layer slips is inversely proportional to the layer thickness as the following 11 :
where h is the layer thickness, μ the shear modulus, ν the Poisson's ratio, φ the angle between slip plane and the layer interface, b the Burgers vector of dislocation, and α is a coefficient representing the extent of dislocation core. Therefore, the larger the layer thickness, the smaller the critical stress is to activate dislocation motion within the layer. The dislocation induced plasticity would increase with the layer thickness, and therefore increase the Green strain and suppress the crack propagation. In SX samples, dislocations can only come from the crack tip due to the high stress concentration, and the number of dislocation sources is limited. In contrast, the available dislocation sources are plenteous in Cu/Nb sample due to the high density of interfacial dislocation networks 9,29 . That is why the Green strain in Cu layer in 20 nm Cu/Nb sample is even larger than that in the SX Cu sample.
Figure 5 shows the atomistic structure for different PX samples at 12 % engineering strain.
Figure 5 (a)-(c) present the samples for PX Cu, Cu layers in 5nm and 20 nm PX Cu/Nb samples, respectively. For PX Cu, the crack propagated within the initial grain firstly. After the tip approach the GBs, the propagation stopped and the crack was blunted at GBs as shown in Figure 5 (a). For the Cu layer in 5 nm PX Cu/Nb sample, after the preexisting crack approached the closest GBs, it continued growing along other GBs. That induced the intergranluar fracture in the thin Cu layer. However, when the layer thickness of Cu layer increased to 20 nm, the intergranular fracture disappeared and crack also was blunted at GBs as shown in Figure 5 (c).
According to the CLS model, the critical stress for dislocation slip within each thin layer in 20nm samples is much lower than that for 5nm samples. Thus, more dislocations can slip at the front of the crack, especially those nucleating from the GBs, to release the elastic energy and blunt the crack tip in the 20nm sample rather than open new surfaces at the front of the crack along the GBs. Figure 6 (c) compares the total Green strain induced in PX samples. The Green strain induced by grain boundaries atoms were not considered on this plot, as the shape, size of grains and the GBs types are exactly the same for all PX samples. It is clear that the Green strain in the PX Cu samples is always the highest one among all cases, followed by the Cu layer in 20 nm Cu/Nb sample. While the single phase PX Nb displays lowest Green strain for the full range of applied strain. This trend is consistent with the crack propagation distance plot in Figure 6 (a).
Since the excess stored elastic energy can either be released by the crack propagation or the dislocation slip, larger Green strains induced by the dislocation slip can suppress the crack propagation and increase the SIF. Thus, the Green strains in Nb layers bonded by the Cu/Nb interfaces are higher than that in the single phase PX Nb, and increasing the layer thickness can facilitate the dislocation activities to generate more Green strains.
In Cu and Nb layer by promoting the intergranular fracture, so the effect of interface may be not be as prominent as that in SX Cu/Nb samples when the layer thickness is small, while the effect of layer thickness on the fracture resistance will not be altered by the GBs. Our findings in this work can provide fundamental understanding of the fracture behavior of MNCs and have implications for the design of nanostructured materials with better fracture resistance.
We perform the MD simulation by using Large-scale Atomic/Molecular Massively Parallel simulator (LAMMPS) codes 30 . Interatomic potential based on Embedded Atom Method Cu-Nb respectively 2,31,32 . These potential have been used widely in the last decade for many studies and provided a lot of insight for understanding the deformation mechanism [33][34][35][36] .
Before experiencing loading, the samples were relaxed by the conjugate gradient method firstly, then equilibrated at 300K for about 40 ps by the Nose/Hoover isobaric-isothermal ensemble (NPT) 27 and the pressure in Z direction was kept zero 37,38 . After relaxation, we load the sample by increasing the stress intensity factor of √ per step based on the fracture mechanics solution for mode-I fracture 21 . During the each loading step, the boundary atoms within 1 nm from the edge in X and Y directions were fixed while other mobile atoms were allowed to relax for 1 ps. The crack tip was recorded to extract the crack length at a specific applied strain. Common neighbor analysis method 32 and Green strain tensor 39 were calculated to characterize the microstructure evolution. Dislocation structures were generated by the dislocation extraction algorithm (DXA) 40 . Finally, the atomistic structures were visualized by the software OVITO 41 .
To determine the Green strain, , for each atom, we firstly calculated the local deformation gradient tensor F for each atom based on the derivative of the relative displacements of the atom's neighbors (1) where X = [X 1 , X 2 , X 3 ] is the original distance vector for atoms to their references before deformation, and x = [x 1 , x 2 , x 3 ] is the new distance vector during deformation. The neighboring atoms should locate within a cutoff radius for 3.5 Angstrom, which can include at least three non-coplanar neighbors for the targeted atom 42 . Then the Green strain is calculated by:
(
where is the identity tensor.