A multiscale FEM-MD coupling method for investigation into atomistic-scale deformation mechanisms of nanocrystalline metals under continuum-scale deformation

Since the potential energy determines forces acting on each atom, the selection of an interatomic potential function significantly affects the results of MD simulations. In nanocrystalline metals, crystal defects, such as stacking faults and deformation twins, are generated during deformation. Therefore, it is critical to select a potential suitable for reproducing these defects. Moreover, for accurate prediction of the properties of metallic materials, taking the local atomic charge into account is required. The many-body potential, which is referred to as the embedded atom method (EAM), is a class of interatomic potential functions based on the density functional theory and the notion of a quasi-atom [47]. For example, the EAM potential proposed by Mishin et al [48] used the ab initio calculations to evaluate the accuracy of parameter fitting for metallic materials with a crystal structure, improving the correspondence between experimental and simulation data. This makes it possible to accurately reproduce basic equilibrium properties, such as the elastic modulus, vacancy formation, and stacking fault energy. The total energy is expressed as the sum of the two-body term that represents the contribution from the two-body potential and the embedding energy:

$$\begin{eqnarray}&&{E}_{\mathrm{tot}}=\displaystyle \frac{1}{2}\displaystyle \sum _{{ij}}^{N}V\left({r}_{{ij}}\right)+\displaystyle \sum _{i}^{N}F\left({\bar{\rho }}_{i}\right),\end{eqnarray} \tag{ 1 }$$

where N is the number of atoms in the system, $V\left({r}_{{ij}}\right)$ denotes the interatomic energy, r_ij represents the interatomic distance between atom i and atom j, $F\left({\bar{\rho }}_{i}\right)$ is the embedded atom energy at the location of atom i with the electron density ${\bar{\rho }}_{i}$ given by

$$\begin{eqnarray}&&{\bar{\rho }}_{i}=\displaystyle \sum _{j\ne i}^{N}\rho \left({r}_{{ij}}\right),\end{eqnarray} \tag{ 2 }$$

where $\rho \left({r}_{{ij}}\right)$ is the atomic density function determined by the interatomic distance r_ij .

The overall effective stress P is given by [49] as,

$$\begin{eqnarray}&&{\boldsymbol{P}}=\displaystyle \frac{1}{V}\left(\displaystyle \sum _{i=1}^{N}{m}_{i}{{\boldsymbol{v}}}_{i}\otimes {{\boldsymbol{v}}}_{i}+\displaystyle \sum _{{\boldsymbol{n}}\in {Z}^{3}}\displaystyle \sum _{i=1}^{N}{{\boldsymbol{r}}}_{i}^{{\boldsymbol{n}}}\otimes {{\boldsymbol{p}}}_{i}^{{\boldsymbol{n}}}\right).\end{eqnarray} \tag{ 3 }$$

The first term on the right-hand side of equation (3) is the contribution from the kinetic energy, which is a function of the velocity v _i and the mass m_i of atom i. The second term on the right-hand side of equation (3) is the Virial tensor with ${{\boldsymbol{r}}}_{i}^{{\boldsymbol{n}}}$ being the coordinate of atoms, ${{\boldsymbol{p}}}_{i}^{{\boldsymbol{n}}}$ the forces by the surrounding atoms, and n ∈ Z³ the three-integer vector representing the shift in the x, y, and z-direction of the periodic systems from the local cell. For example, ${{\boldsymbol{r}}}_{i}^{{\boldsymbol{n}}}\,+\,{\boldsymbol{Ln}}$ with L being the periodic cell vector means the position of atom i is shifted by ${\boldsymbol{Ln}}$ from that in the local cell. The restriction on the summation of the virial tensors is imposed so that only one periodic image of interactions of atoms is considered. With this consideration, only one self-interaction with a mirror atom, which is the same atom generated by periodic boundaries, will be taken into account for a computation.

For a domain V with surface A in the absence of volumetric contribution and neglecting inertia contribution, the principle of virtual work is written using the Voigt notation, i.e., collecting components of a tensor in a column, as

$$\begin{eqnarray}&&{\iiint }_{V}\delta {{\boldsymbol{\varepsilon }}}^{T}{\boldsymbol{\sigma }}{dV}-{\iint }_{A}\delta {{\boldsymbol{u}}}^{T}{\boldsymbol{t}}{dA}=0,\end{eqnarray} \tag{ 4 }$$

where δ ε is the virtual strain, σ is the stress, δ u is the virtual displacement, and t is the surface traction. The interpolation in terms of the columns of the nodal displacements a and the virtual nodal displacements δ a enables expressing the actual and virtual displacements u and δ u and the actual and virtual strains ε and δ ε as

$$\begin{eqnarray}&&{\boldsymbol{u}}={\boldsymbol{N}}{\boldsymbol{a}},\quad {\boldsymbol{\varepsilon }}={\boldsymbol{B}}{\boldsymbol{a}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\delta {\boldsymbol{u}}={\boldsymbol{N}}\delta {\boldsymbol{a}},\quad \delta {\boldsymbol{\varepsilon }}={\boldsymbol{B}}\delta {\boldsymbol{a}},\end{eqnarray} \tag{ 6 }$$

where N is the column of the finite element shape functions and B is the matrix containing the spatial derivatives of the shape function. Substituting equation (6) and taking into account that equation (4) has to hold for all virtual displacements, ∀ δ a , lead to

$$\begin{eqnarray}&&{\iiint }_{V}{{\boldsymbol{B}}}^{T}{\boldsymbol{\sigma }}{dV}-{\iint }_{A}{{\boldsymbol{N}}}^{T}{\boldsymbol{t}}{dA}=0.\end{eqnarray} \tag{ 7 }$$

The linear elastic constitutive equation for the stress is represented in an incremental form as,

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}(k)={\boldsymbol{D}}\left({\boldsymbol{\varepsilon }}(k)-{\boldsymbol{\varepsilon }}(k-1)\right)+{\boldsymbol{\sigma }}(k-1),\end{eqnarray} \tag{ 8 }$$

where D is the tensor that provides the relationship between the strain and the stress, ε (k) and σ (k) are the strain and the stress at step k, respectively. Substituting equation (8) into equation (7) gives the discretized equilibrium equation:

$$\begin{eqnarray}&&\left({\iiint }_{V}{{\boldsymbol{B}}}^{T}{\boldsymbol{D}}{\boldsymbol{B}}{dV}\right){\boldsymbol{a}}(k)={\iiint }_{V}{{\boldsymbol{B}}}^{T}{\boldsymbol{D}}{\boldsymbol{\varepsilon }}(k-1){dV}-{\iiint }_{V}{{\boldsymbol{B}}}^{T}{\boldsymbol{\sigma }}(k-1){dV}+{\iint }_{A}{{\boldsymbol{N}}}^{T}{\boldsymbol{t}}{dA}.\end{eqnarray} \tag{ 9 }$$

The integrals in equation (9) are usually computed using a numerical integration scheme, e.g., the Gauss quadrature. In the Gauss quadrature, an integral is evaluated using an appropriate number of specifically positioned IPs with the corresponding weights. Finally, the global discrete equilibrium equation equation (9) is expressed as

$$\begin{eqnarray}&&{\boldsymbol{K}}{\boldsymbol{a}}(k)={\boldsymbol{f}}(k-1),\end{eqnarray} \tag{ 10 }$$

where K is the stiffness matrix and f is the force acting on the domain. Nodal displacements a are solved using equation (10) with a given K and boundary conditions.

In our FEM-MD multiscale simulations, we utilize the notion of CBR, a hierarchical hypothesis for crystalline solids, to relate the positions of atoms with a continuum model deformation. CBR states that the positions of the atoms within the crystal lattice follow the macroscopic strain of the corresponding continuum medium [50]. Let the deformation gradient tensor of a point of a continuum model be F . Then, the infinitesimal elemental vector d X in the reference configuration and the corresponding vector d x in the current configuration are related as follows,

$$\begin{eqnarray}&&d{\boldsymbol{x}}={\boldsymbol{F}}d{\boldsymbol{X}}.\end{eqnarray} \tag{ 11 }$$

A crystalline solid's lattice vector is considered an elemental vector. We chose, here, the classical CBR as the considered deformation is not very large and the size of the MD cells is small compared to the spatial variation of macroscale deformation; thus, the discrepancy due to the low order is negligible. Another reason is that classical CBR suffices for the problem of crystalline bulk as we do not consider a localized deformation on MD unit cells. Figure 1 schematically illustrates the coupling of the deformation of a continuum model and a representative atomistic cell using CBR in the case of shear deformation.

Zoom In Zoom Out Reset image size

The method of coupling the MD model with the FEM using the concept of CBR is explained in this section. Figure 2 represents the multiscale computational algorithm employed in this study. Firstly, we prepare a representative MD simulation cell whose shape is described by a second-order shape tensor H . Each computational cell used in MD is virtually related to a Gauss quadrature IP of a finite element model, as shown in figure 3. Each MD simulation representing the corresponding integration point of the finite element is independent of each other. Direct communication between the MD cells is not considered as each cell is made of a locally periodic nano-scale structure that has a large distance from other cells at a macroscale. The atom position vector r _i (k) for an atom i in a representative unit cell at time increment k can be expressed as

$$\begin{eqnarray}&&{{\boldsymbol{r}}}_{i}^{\mathrm{IP}}(k)={{\boldsymbol{H}}}^{\mathrm{IP}}(k){{\boldsymbol{r}}}_{i}^{\mathrm{IP}}(0)\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\mathrm{with}\quad {{\boldsymbol{H}}}^{\mathrm{IP}}(k)={{\boldsymbol{F}}}^{\mathrm{IP}}(k){{\boldsymbol{H}}}^{\mathrm{IP}}(0),\end{eqnarray} \tag{ 13 }$$

where F ^IP(k) is the deformation gradient tensor of IP at time increment k in the finite element model and H ^IP(0) is the initial shape of the representative MD unit cell. The shape tensor H ^IP(k) prescribes the shape of the unit cell at step k. However, F ^IP(k) is, in general, not symmetric and has nine independent components. The rotation component of the deformation gradient tensor needs to be removed so as not to allow a rigid body rotation of MD unit cells. To this end, we introduce the QR decomposition to eliminate the rotation components. Using the QR decomposition by the Gram-Schmidt orthogonalization, the deformation gradient tensor F ^IP(k) is decomposed into the orthogonal rotation tensor Q ^IP(k) and the upper triangular stretch tensor R ^IP(k):

$$\begin{eqnarray}&&{{\boldsymbol{F}}}^{\mathrm{IP}}(k)={{\boldsymbol{Q}}}^{\mathrm{IP}}(k){{\boldsymbol{R}}}^{\mathrm{IP}}(k).\end{eqnarray} \tag{ 14 }$$

The upper triangular stretch tensor R ^IP(k) satisfies the periodic boundary conditions on the orthogonal basis of the unit cell. Therefore, instead of using the full F ^IP(k), we use R ^IP(k) to apply the deformation to the unit cell:

$$\begin{eqnarray}&&{{\boldsymbol{H}}}^{\mathrm{IP}}(k)={{\boldsymbol{R}}}^{\mathrm{IP}}(k){{\boldsymbol{H}}}^{\mathrm{IP}}(0).\end{eqnarray} \tag{ 15 }$$

In an incremental form, equation (15) can be rewritten as

$$\begin{eqnarray}&&{{\boldsymbol{H}}}^{\mathrm{IP}}(k)={\rm{\Delta }}{{\boldsymbol{R}}}^{\mathrm{IP}}(k){{\boldsymbol{H}}}^{\mathrm{IP}}(k-1),\end{eqnarray} \tag{ 16 }$$

where the increment Δ R ^IP(k) is obtained using the values at the previous step k − 1 and the current step k as,

$$\begin{eqnarray}&&{\rm{\Delta }}{{\boldsymbol{R}}}^{\mathrm{IP}}(k)={{\boldsymbol{R}}}^{\mathrm{IP}}(k){\left({{\boldsymbol{R}}}^{\mathrm{IP}}(k-1)\right)}^{-1}.\end{eqnarray} \tag{ 17 }$$

Therefore, new atomic coordinates are estimated based on the previous atom position vector ${{\boldsymbol{r}}}_{i}^{{IP}}(k-1)$:

$$\begin{eqnarray}&&{{\boldsymbol{r}}}_{i}^{\mathrm{IP}}(k)={\rm{\Delta }}{{\boldsymbol{R}}}^{\mathrm{IP}}(k){{\boldsymbol{r}}}_{i}^{{IP}}(k-1)\end{eqnarray} \tag{ 18 }$$

After equation (18) is prescribed, the coordinates of the atoms within the unit cell are determined by minimizing the total potential energy to approximate the equilibrium state. Next, the information about the atom coordinates and forces is used for the evaluation of the stress and the analysis of the crystal structure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At the continuum macroscale, a cubic model with an edge length of 0.4 mm is considered. This macroscopic domain is discretized with a single isoparametric linear eight-node hexahedral element with eight IPs. The macroscopic element with the nodes and the IPs are depicted in figure 4(a).

Zoom In Zoom Out Reset image size

Two non-uniform tensile loading cases, which are described in the following parts, are applied to the macroscopic domain. In both loading cases, three faces on the left, the bottom, and the front were fixed with sliding allowed; that is, node 1 was fully fixed, node 2 was fixed in the y- and z-directions, node 5 in the x- and y-directions, node 3 in the x- and z-directions, node 7 in the x-direction, node 6 in the y-direction, and node 4 in the z-direction. In loading case 1 (figure 4(b)), we prescribe different values of displacements to nodes 2, 4, 6, and 8 on the right side of the element, shown in figure 4(b). Displacements in the x-direction of 0.05 mm, 0.06 mm, 0.07 mm, and 0.08 mm were given to nodes 2, 4, 6, and 8, respectively. In loading case 2 (figure 4(c)), only nodes 2 and 4 were given displacements of 0.05 mm in the x-direction. For both loading cases, the total displacements were applied in 30 steps. After the simulation of one step in FEM was completed, the deformation state of all nanocrystalline models corresponding to different IPs was updated accordingly. The prescribed boundary conditions on the nodes are summarized in table 1.

At the MD scale, a nanocrystalline aluminum model is considered, which consists of four hexagonal grains with almost identical grain sizes; see figure 5. The MD model was generated using Atomsk [51]. The unit cell has a cubic shape, with each side being approximately 30 nm; see appendix A on how the size of the unit cell was chosen. We consider that the columnar-structured nanocrystalline model is sufficient to demonstrate the FEM-MD method, at this moment, not aiming at an in-depth study of the detailed deformation mechanisms of nanocrystalline aluminum. Nevertheless, grain size has been selected such that the typical behavior of nanocrystalline aluminum is well represented; see appendix A. The crystal orientation of the grain labeled (a) is defined as x [−1 1 0] y [−1−1 2] z [1 1 1]. The grains labeled (b), (c), and (d) have the same crystal orientation in the z-direction but different x- and y-directions obtained by the counterclockwise rotation around the z-axis by 30°, 60°, and 90°, respectively, from the crystal orientation of grain (a); see table 2. This configuration is denoted as nanocrystalline model A. The distinct crystal orientations generate the grain boundaries appearing as gray-colored regions in figure 5. Another nanocrystalline model B has a different crystal orientation pattern. In model B, the crystal orientation of each grain was rotated 45 degrees around the z-axis from the original orientation in nanocrystalline model A; see table 2.

Zoom In Zoom Out Reset image size

The MD simulations were performed in LAMMPS [43]. Periodic boundary conditions were applied to all MD cell boundaries. Mishin's EAM was chosen as the interatomic potential [48]. After the generation of the nanocrystalline models, initial relaxation of the unit cells was performed. It took 20 000 steps (20ps) for each computation in the NPT ensemble, in which the pressure 0 Pa and the temperature 300 K are kept constant. The unit cell of the nanocrystalline model was assumed to be embedded in each IP. The deformation information was passed to the computational cell of the MD method through the connection between each IP and the MD model. For simplicity, the position vector defined in equation (18) was applied only in the x-direction to the MD unit cells. The y- and z-directions were pressure-controlled. The MD computations to deform each computational cell were conducted at a constant temperature of 300 K and a constant pressure of 0 Pa in the y- and z-directions. Each computation of a deformation phase took 500 steps, which is 0.5 ps in the time unit. After the unit cell had finished the prescribed deformation process, a relaxation phase was inserted to make the systems close to equilibrium states. The relaxation of the system was performed for 35,000 steps (35 ps) at constant pressure with the temperature set at 300 K; see appendix B for the study on the length of relaxation time. After the computation of the relaxation had finished, we conducted the Common Neighbor Analysis (CNA) to observe the nucleation and the growth of defects such as deformation twins and stacking fault [52]. These processes were repeated until all the computation steps prescribed on the FEM side were done.

The macroscopic deformation of loading case 1 is shown in figure 4(b). Four distinct deformation gradient values at IPs 1 and 2, 3 and 4, 5 and 6, and 7 and 8 were obtained for this loading. The effective stress-strain responses resulting from the MD simulations are shown in figure 6. It is argued that strain softening occurred at all the IPs after reaching the yield strain of approximately 0.025. The trend of the yield strains agrees well with that obtained in the experimental tensile test of nanocrystalline aluminum [53], both of which have a yield strain between 0.02 and 0.025. The yield stress differs between the MD simulation and the experiment due to the high strain rate in the MD simulation.

Zoom In Zoom Out Reset image size

We also assessed the ratio of hexagonal close-packed (hcp) atoms with respect to the total number of atoms in the model at every FEM step after the relaxation phase in MD computations to investigate the nucleation of crystal defects. The number of hcp atoms is considered in the discussion primarily to demonstrate that the stress-strain response correlates with the onset and evolution of stacking faults or twins. As shown in figure 7, an increasing trend in the number of hcp atoms is observed after passing the yield strain value. The results also show that the hcp ratio at IPs 1 and 2 was, on average, higher than at the other IPs. The highest hcp ratio at IPs1 and 2 was 0.0190, whereas the lowest value was 0.0118 at IPs7 and 8. Figure 8 shows the crystal structure in cross-section at z = 20 mm at selected FEM steps after relaxation. Multiple deformation twins crossing the central grain in the diagonal direction were observed in the nanocrystalline MD models at IPs 1 and 2, as shown in figures 8(a) and (b). At IPs 3 and 4, on the other hand, only a few crystal defects in the nanocrystalline MD models occurred, which were also confirmed by a low hcp ratio in the plastic deformation stage. Here, we observed that the stacking faults that existed at the applied strain ε_xx = 0.053 in figure 8(c) vanished at ε_xx = 0.053 in figure 8(d). While some stacking faults disappeared, other defects were nucleated at the same time. This phenomenon also occurred for the stacking faults and deformation twins at other locations at IPs 7 and 8, as seen in figure 8(h).d Further investigation into the behavior of defects was conducted by studying the crystal structure changes during the relaxation, shown in figures 9 and 10. These snapshots exhibit the appearance, annihilation, and migration of crystal defects in the course of the relaxation. At the beginning of the relaxation figure 9(a), no crystal defects are present, whereas the defects represented by red lines are confirmed in figure 9(d). Similar observations are made in figures 9(e)–(h). The snapshots in figures 9(e)–(h) exhibit the formation of deformation twins at the early stage of the relaxation, and the size of the twins grows in the course of the relaxation. The crystal structure change during the next relaxation phase following the next loading step is shown in figure 10. The annihilation of defects is observed in figures 10(b)–(d). Furthermore, detwinning is observed in figures 10(f)–(h). The onset of various crystal structure changes, those shown in the snapshots in figures 9 and 10, typically originates from the emission of defects from the grain boundaries as indicated in figure 11, see also the supplementary video. Compared to the initial configuration shown in figure 5, the grain boundaries after loading are not neatly aligned and are somewhat jagged. Previous studies have already revealed that these disordered grain boundaries, called nonequilibrium grain boundaries (NEGBs), induce atomic shuffling and rearrangement in the MD simulation [54, 55]. Similar observations have also been confirmed in the experimental investigation into nanocrystalline materials as well [56, 57]. The present work considers more complicated grain structures and deformation history than those in the previous studies and thus suggests that the interaction between NEGBs could result in the onset and annihilation of deformation twins during relaxation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

From the deformation shown in figure 4(c), four distinct deformation patterns were obtained at IPs 1 and 2, 3 and 4, 5 and 6, and 7 and 8, as well as in case 1. The effective MD stress-strain curves for the xx component obtained at corresponding IPs are shown in figure 12. While strain softening was confirmed at IPs 1-4, the stress response remained elastic at IPs 5-8 because of the small deformation prescribed on the unit cells. The result of figure 13 indicates that the number of hcp atoms increased at IPs 1-4, which is attributed to the onset or growth of defects after reaching the yield point. The ratios of hcp atoms for IPs 5-8 are not plotted since they remained almost zero during the whole simulation. The snapshots of crystal structure in figure 14 show that IPs 1 and 2 contained deformation twinning under deformation. On the other hand, no deformation twin was found at IPs 3 and 4, although stacking faults and some other defects were observed in figure 14(c). It is also noteworthy that in both configurations, the grain boundaries are no more as straight as they are in the initial configuration. As discussed in loading case 1, it is suggested that NEGBs also came into play in the formulation and annihilation of stacking faults and deformation twins under this loading pattern.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The deformation experienced at each IP is identical to loading case 1. The effective stress-strain curves shown in figure 15 and the result of CNA provided in figure 16 showed the different responses from those obtained from the simulation for nanocrystalline model A. Figure 15 confirmed that the maximum stress was reached at a strain of approximately 0.04, similar to the other two cases. Notably, the increase in the ratio of hpc atoms in figure 16 was smaller than in the case of the nanocrystalline model A as in figure 7. In this case, the hcp ratio at a strain of 0.05 is below 0.005. On the other hand, looking at figures 7 and 13, one can see that in the other two cases, the hcp ratio at a strain of 0.05 is higher than even 0.01. Since the increase in the number of hcp atoms is associated with the evolution of crystal defects, this observation indicates orientation-dependent defect activities, which have been observed in experiments and simulations on conventional metals and nanocrystalline metals. We confirmed that more crystal defects were generated in nanocrystalline model A in figure 8 compared to figure 17. However, onsets of stacking faults were still observed in nanocrystalline model B. No deformation twin was found in figures 17(a)–(d), which can also be explained by the orientation dependency. Furthermore, NEGBs are also seen in figure 17 in comparison with the initial configuration in figure 5, which infers the interplay of NEGBs in the case of nanocrystalline model B as well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To summarize, the developed FEM-MD multiscale method allows observation of the evolution of various dislocation patterns, including stacking faults and deformation twins, in the atomistic nanocrystalline model as a response to the deformation history originating from the macro-scale FEM. Regarding validation, the yield strain and crystal structure changes show good qualitative agreement with experimental observations. The variation of the results among different IPs is the most pronounced when some IPs do not show any defect evolution while others do. Furthermore, the extent of generation of stacking faults and deformation twins varies between IPs, exhibiting different dislocation nucleation patterns depending on the local macroscale deformation scenarios. Another remark is that different realistic microstructural evolution patterns are seen depending on the combination of different deformation scenarios and crystal orientations, which provides guidelines for the optimum material design of nanocrystalline materials. On the other hand, one of the limitations lies in that certain extreme loading conditions, such as twisting and bending, cannot be modeled in the MD simulation using the present implementation. Another limitation is the computational cost if a large model is employed. One of the solutions to this problem is to utilize k-means clustering on integration points with similar deformation patterns, leading to a reduction in computational cost. The employment of k-means clustering has already been demonstrated in some studies on the homogenization method [58–60]. Nevertheless, the present method can evaluate the dislocation dynamics of nanocrystalline models that contain multiple grains based on realistic macroscale deformation history.