Molecular dynamics simulations of mechanical behavior in nanoscale ceramic–metallic multilayer composites

ABSTRACT The mechanical behavior of nanoscale ceramic–metallic (NbC/Nb) multilayer composites with different thickness ratios is investigated using molecular dynamics (MD) simulations. Based on the obtained stress–strain behavior and its dependence on temperature, strain rate, and loading path, the flow stress for the onset of plasticity is identified and modeled based on the nucleation theory, and the in-plane yield loci for different layer thicknesses are constructed. The results are used to establish the plastic flow potential for developing a continuum viscoplastic constitutive model for potential use in large-scale applications. GRAPHICAL ABSTRACT IMPACT STATEMENT Using MD simulations, we provide new understandings of the mechanical behavior of ceramic–metallic nanolaminates by constructing the yield loci and proposing a plastic flow potential under parallel-to-interface biaxial loading conditions.


Introduction
High strength and lightweight materials are the center of attention due to their extensive applications in various industries such as aerospace, automotive, and oil. Ceramic-metallic nanolaminates (CMNs) are potential candidates for such applications as they present high strength and ductility, high wear and corrosion resistance, and the ability to maintain superior properties at high temperatures [1,2]. The optimized design of CMNs is achieved by selecting the adequate combination of individual layer thickness, thickness ratio of the CONTACT Mohsen Damadam mohsen.damadam@wsu.edu Spokane and college St, Voiland college of engineering and architecture, ETRL 236, Pullman, WA 99163, USA multilayer constituents, interface type and structure, and mechanical properties of the individual layers [2][3][4]. The properties of the interface largely influence the plastic deformation of CMNs when the thickness of the individual layers is in the order of nanometers [1,3,5]. Interfaces are categorized into two groups according to their structure: coherent and semi-coherent interfaces. The interface structure, in turn, is dependent on the adjoining layers' crystal structure and orientation. In coherent (transparent) interfaces, atoms on both sides have strictly one-to-one correspondence and registry, such that dislocations in one layer can pass to the neighboring layer by overcoming the coherency stresses at the interface [6,7]. Coherent interfaces commonly exist between layers with identical lattice structures, very small lattice mismatch, and very small layer thickness. The semicoherent interface comprises coherent regions separated by intrinsic interface dislocations. It is a superposition of an intrinsic interface dislocation network with a perfectly coherent interface, where misfit dislocations relax the far-field coherency stress of the coherent interface. Localized and extended nodes can form at the intersection of misfit dislocations in the interface [8][9][10][11]. Semicoherent interfaces can exist between layers that either share the same lattice structures and orientation, such as Cu-Ni and Cu-Ag [12], or have completely different crystal structures, such as Cu-Nb and Mg-Nb [2,13]. Some semi-coherent interfaces, such as those between dissimilar crystal types or the same crystal type but with very large lattice mismatch, have small coherent regions and very densely arranged interface dislocations. The strong overlapping of dislocation cores made the structural features of the interface difficult to identify [9,10]; these interfaces have, therefore, been traditionally classified as 'incoherent (opaque) interfaces'.
Apart from the interface, both the thickness of the individual layers and the thickness ratio between the two layers have a considerable influence on the mechanical response of CMNs which has been investigated experimentally under various loading conditions including compression [14][15][16], tension [17], and nanoindentation [18][19][20]. The results showed a reduction in strength and hardness upon increasing the layer thickness. Furthermore, plastic co-deformation of metal and ceramic in CMNs is possible when both the absolute thickness of the ceramic layer and the thickness ratio between the ceramic layer and the metal layer are small [4,[21][22][23]. Under uniform compression, and in the absence of propagating dislocations inside the individual layers of NbC/Nb multilayers, dislocations nucleate from localized nodes where the local strain is the highest associated with the core of the dislocation intersection. In addition, if NbC/Nb nanolaminates are under nanoindentation, depending on the location of the indenter, new dislocations nucleate either from the localized nodes or the misfit dislocation lines in the interface [8].
Previously, specific yielding and subsequence plastic deformation mechanisms of CMNs have been investigated in detail. However, effective and efficient prediction of the strength of a bulk CMN still demands a more thorough understanding of the yielding behavior of the CMN under a more generalized framework. In this study, we construct complete yield loci of NbC/Nb nanolaminates with varying individual layer thickness/thickness ratios. The plastic deformation of NbC/Nb nanolaminates was investigated under uniaxial tension and compression using atomistic simulations. The yield loci were fitted with a function representing a plastic flow potential for NbC/Nb nanolaminates. The plastic flow potential of CMNs can be combined with a continuum viscoplastic framework and integrated to a finite element code for large-scale structural analysis to be able to study the thermo-mechanical behavior of bulk CMNs.

Simulation set-up and results
Molecular dynamics simulations of NbC/Nb nanolaminate composites are carried out using LAMMPS [24] under various uniaxial and in-plane biaxial loading conditions. First, simulations were carried out for a range of temperatures (10-1000 K) and strain rates (10 7 -10 9 s −1 ) under uniaxial tension with the purpose of establishing the flow stress and its dependence on temperature and strain rate. It is observed that at 10 K, the nucleation stress dependence on the strain rate is almost negligible. Furthermore, at this low temperature thermal fluctuations are negligible, which, in turn, provide conditions for a more clear visualization of the deformation mechanisms. Next, and in order to establish the shape of the yield surface, a series of biaxial loading conditions were performed at constant temperature (10 K) and constant strain rate (3 × 10 8 s −1 ) for different layer thicknesses. From these sets of simulations, we establish a relation describing the yield surface and its dependence on temperature and strain rate.
Throughout all the simulations, a second nearestneighbor modified embedded-atom method (2NN-MEAM) [25,26] interatomic potential was implemented for both Nb and NbC layers. This interatomic potential is capable of reproducing different physical, elastic, thermal, and interfacial properties for NbC and Nb that are in proper agreement with experimental information [8,27]. Periodic boundary conditions were applied in all directions to study the behavior of the bulk material.
The difference between the crystallographic structures and lattice constants of Nb and NbC, that is, BCC and rock salt (B1), respectively, results in the formation of lattice misfit dislocations at the NbC/Nb interface. The atomistic structure involves a single layer of each ceramic (NbC) and metal (Nb) brought together and rotated to get the desired crystallographic orientation based on the Baker-Nutting orientation relationship between the two lattices, that is, {001} Nb ||{001} NbC and 110 Nb || 110 NbC . Additional strains were applied to both layers, creating a strainednatural dichromatic pattern (strained-NDP), to assure the periodicity in lateral directions [13,28]. The geometry of the NbC/Nb bilayer structure is illustrated in Figure 1. The interface is parallel to the x-z plane, while the y axis is perpendicular to the interface.
The strains applied along the x and z directions result in a structure with both layers having the same sizes in the x and z directions. The theory of elasticity is used to find the dimensions of the structure in the x and z directions [13,28]. Assuming that the stress normal to the interface (x-z plane) is zero, the minimum applied strains in the x and z directions, for a defined individual layer thickness, are obtained using Hooke's law and the equilibrium equations, while the shear stress at the interface is assumed to be negligible: 2C NaCl 12 ε NaCl where C ij 's are the components of the stiffness tensor, and the thickness of NbC and Nb is defined as t NaCl y and t bcc y , respectively. N NaCl,bcc x,z and α NaCl,bcc x,z are, respectively, the number of atomic planes and lattice spacing of the NbC or Nb layers along the x or z direction. The minimum strains applied along the x and z directions for different thickness ratios are found to vary between 10 −5 and 10 −4 such that the NbC layer is stretched and the Nb layer is compressed to obtain the registry between the layers.
Equations (1)-(4) reveal that for samples with different thicknesses but the same thickness ratio, the same amount of strain is needed to be applied to the layers. The dimensions of the structure are determined to be 20.65 nm in the x and z directions, corresponding to 44 unit cells for NbC and 43 unit cells for Nb. The number of atoms in the simulation box was between 242,000 and 285,000 depending on the individual thickness of each layer. For all the temperature cases studied, the initial configuration was equilibrated, first by applying energy minimization at 0 K until the maximum force on any atom was smaller than 1 × 10 −10 N. Then dynamic relaxation was applied for 20 ps through NPT ensemble at zero pressure. In order to verify that the relaxation time was enough, simulations were also performed up to 40 ps and it was observed that similar interfacial misfit dislocation networks were obtained. The relaxed structure was then subjected to a uniaxial tensile or compressive loading parallel to the interface plane with a constant strain rate.
Three cases of constant strain rates were performed, 3 × 10 7 , 3 × 10 8 , and 10 9 s −1 , which are typical values for MD simulations [29,30]. The stress-strain responses for the CMNs subjected to uniaxial loading were obtained and the yield points were identified. Figure 2(a,b) shows typical results for the stress-strain curves, for the case of 10 K and strain rate of 3 × 10 8 s −1 , for NbC/Nb multilayers with a total bilayer thickness of 10 nm and various ceramic-metallic thickness ratios, subjected to tensile and compressive loadings in the x direction, respectively. Increasing the metal thickness results in a decrease in the elastic modulus, leading to a significant decrease in the yield strength. Under uniaxial tension with a thickness ratio of one (5 nm NbC/5 nm Nb, red solid curve in Figure 2(a)), dislocations nucleate and propagate ( Figure  2(c)) in the Nb layer without any evidence of ductility in the NbC layer. Cracks form in the NbC layer shortly after the yield of the Nb layer, and the CMN fails (marked by 'x' in Figure 2(a)). However, ductile behavior is observed in the NbC layer by decreasing the ceramic layer thickness, evidenced by the presence of two peak stresses (green short dash curve and blue long dash curve in Figure 2(a)) and the observation of dislocations propagation in the NbC layer (Figure 2(d)). Under uniaxial compression (Figure 2(b)), dislocations nucleate in the Nb layer and transmit to the NbC layer without any crack initiation. The characteristics of the dislocations within the NbC/Nb interface such as dislocation type and their Burgers vectors were comprehensively studied by Salehinia et al. using atomistically informed Frank-Bilby theory [5]. It is noted that the snapshots in Figure 2(c,d) are captured in OVITO [31].
It is noted that the yield strengths predicted in this work are for the cases of high strain rates (10 7 -10 9 s −1 ), and thus are in the order of 10-20 GPa, which far surpasses the typical range of experimental observations for much lower strain rates [32]. This is due to the limitation in the timescale of the MD simulations. However, the MD results can be used to formulate, based on the nucleation theory, a MD-based constitute equation that can predict the flow stress for a wide range of strain rates and temperatures. According to classical nucleation rate theory [33], the rate of dislocation nucleation from the NbC/Nb interface is written as where r* is the critical radius of the dislocation nucleus, b is the Burgers vector, ν is the Debye frequency, n 0 is the total number of nucleation sites, Q* andˆ are the activation energy in the absence of applied stress and activation volume of the nucleation event, and σ is the resolved shear stress. Therefore, a change in the strain rate (and therefore the nucleation rate) by a factor of 10 11 (from 0.001 s −1 typical for experiments [16] to 3 × 10 8 s −1 ) inevitably leads to a significant change in σ . This change can be estimated based on Equation (5) in terms of strain rate and temperature [34]. Based on the nucleation theory, the following constitutive equation for flow stress in nano-layers was derived in [6].
whereσ is a function of the strain rate, layer thickness (here Nb), and a linear function of temperature. The activation parameters can be found by fitting Equation (6) with stress vs temperature results from MD simulations as shown in Figure 3(a) (αβ = 1, v D = 1.3e13 s −1 , h = 5 nm, S = 1), in which the activation volume and nucleation barrier are calculated as: where b is the Nb Burgers vector. Using the nucleation stress model, the predicted flow stresses are comparable to the MD results as can be deduced from Figure  3(b). Thus, this MD-based model can be used to correlate MD results with experimental data for low strain rates. However, one would still expect the stresses to be higher when compared to experiments since the simulations did not include initial defects inside the lattice such as pre-existing dislocations [5].
In order to apply biaxial loading, first, the structure is loaded in the x direction, up to a stress below the yield stress in that direction by applying a constant strain rate of 3 × 10 8 s −1 . Then, the stress is maintained constant (zero strain rate) along the x direction and loading is then applied along the z direction by using the same constant strain rate 3 × 10 8 s −1 . Then, the corresponding stress-strain curves are obtained, from which the yield stress can be determined. During the entire loading process, the normal stress in the y direction is maintained at zero. Typical stress-strain curves for such loadings on two structures, that is, 5 nm NbC/5 nm Nb (red dashed lines) and 3 nm NbC/7 nm Nb (green solid lines), are presented in Figure 4(a,b) for stress in the x and z directions, respectively. It can be seen that by applying the tensile stress up to 12 GPa and compressive stress up to −14 GPa in the x direction and then applying the tensile and compressive stress in the z direction, yielding occurs at 16.5 and −16 GPa, respectively. Hence, the pairs (12, 16.5) and (−14, −16) constitute two distinct points on the yield locus (shown in Figure 4(c)). The red and green curves in Figure 4(a) undergo a decrease and increase, respectively, in strain along the x direction after starting the load in the z direction, which is due to the effect of Poisson's ratio. The complete yield locus was constructed following the same loading path, where different in-plane biaxial loadings such as tension-tension, tension-compression, compression-compression, and compression-tension were applied on the NbC/Nb nanolaminates. Figure 4(c) presents MD-developed convex yield loci obtained by maintaining a constant total bilayer thickness and changing the thickness ratio by decreasing the NbC layer thickness. It is observed that decreasing the NbC layer thickness results in decreasing the yield stress and tension/compression asymmetry. Furthermore, the yield stress when loading uniaxially in the x direction is equal to the flow stress when loading in the z direction for each yield locus. This is attributed to the cubic symmetry along the x and z directions that leads to a symmetric interfacial misfit dislocation pattern. Figure 4(d) presents the MD-developed yield loci of NbC/Nb multilayers where the thickness ratio was 1 and the bilayer thickness was varied from 10 to 20 nm. It is noted that for these two cases, increasing the total thickness of the structure has a slight effect on the overall size of the yield locus, which implies that the critical dislocation nucleation stress is thickness independent for a constant thickness ratio, unless the thickness is extremely small (1-2 nm), where the interaction between interface dislocations on adjacent interfaces becomes significant. When comparing the case of '10 nm NbC/10 nm Nb' shown in Figure 4(d) to the case of '1 nm NbC/9 nm Nb' shown in Figure 4(c), one can conclude that the thickness of the NbC affects the size of the yield locus; increasing the thickness of the NbC layer, while maintaining the thickness of the Nb layer constant, increases the size of the yield locus.

Development of the plastic flow potential
The classical theory of plasticity assumes the existence of a plastic flow potential function that describes the plastic strain rate tensor from an associated flow rule. At the macro scale, the behavior of larger CMN structures can be described using a viscoplastic framework that can be developed utilizing the knowledge gained from the MD simulations. Based on the plasticity theory, yield function is a proper replacement for plastic flow potential [35]. Hence, a general function of a yield locus that describes the anomalous behavior and accounts for the kinematic hardening is introduced to describe the initiation of the plastic deformation in the NbC/Nb nanolaminates [6,36]. Here, we adopt a similar form for the plastic flow potential and incorporate the results obtained in the above section for the flow stress in Equation (8), yielding where σ 11 and σ 22 are the principal stresses,σ is given by Equation (6), and C, α 1 , σ * 1 , α 2 , α * 2 , α 3 , m, and β are the parameters to be determined by fitting the yield function on a particular MD-developed yield locus ( Figure  5(a-c)).  In Equation (8), the term βσ determines the size of the yield locus, while the other terms describe its shape and position in the stress space. The parameters σ * 1 and σ * 2 act as back stress parameters controlling the position of the yield locus [35,36]. Table 1 lists these parameters for the NbC/Nb samples with the same bilayer thickness and different thickness ratios. All fitting parameters are maintained unchanged, except for σ * 1 , σ * 2 , and β. According to the results shown in Table 1, the parameters C, α 1 , α 2 , α 3 , and m, which determine the shape of the yield surface, are constants and independent of the layer thicknesses, while the parameters σ * 1 , σ * 2 , and β are thickness dependent. Decreasing the NbC thickness leads to an increase in the internal stresses (σ * 1 , σ * 2 )) and a decrease in β, which in turn describes the translation Figure 6. MD-based yield loci proposed for NbC/Nb nanolaminates at different thickness ratios. and shrinkage of the yield locus, Figure 5(a-c) and Figure  6. It can be seen from Figure 5 that the proposed function provides a good fit to the MD results. Furthermore, the β values, which are obtained according to Equation (6) (for 5, 7, and 9 nm Nb layer thicknesses) and by fitting Equation (8) with the MD results, suggest that the size of the yield locus is more affected by the NbC layer thickness than by the thickness of the Nb layer.

Conclusions
Molecular dynamics simulations were performed on nanoscale ceramic-metallic multilayer composites under uniaxial and biaxial in-plane loading conditions. Yield loci were obtained for NbC/Nb multilayers at different thickness ratios. The results show that by increasing the Nb layer thickness, the yield strength of the system reduces and consequently leads to the shrinkage of the yield locus. It has also been shown that increasing the Nb layer thickness results in lower tension/compression asymmetry in the yield loci of NbC/Nb nanolaminates. Finally, a function for the plastic flow rule was obtained from the MD-developed yield loci. This flow potential can be further implemented in viscoplastic models to bridge to larger size scales.

Disclosure statement
No potential conflict of interest was reported by the authors.