Spatial distribution of local elastic moduli in nanocrystalline metals

Elastoplastic properties of nanocrystalline metals are non-uniform on the scale of the grain size, and this non-uniformity affects macroscopic quantities as, in these systems, a significant part of the material is at or adjacent to a grain boundary. We use molecular dynamics simulations to study the spatial distributions of local elastic moduli in nano-grained pure metals and analyze their dependence on grain size. Calculations are performed for copper and tantalum with grain sizes ranging from 5-20nm. Shear modulus distributions for grain and grain-boundary atoms were calculated. It is shown that the non-crystalline grain boundary has a wide shear-modulus distribution, which is grain-size independent, while grains have a peaked distribution, which becomes sharper with increasing grain size. Average elastic moduli of the bulk, grains, and grain boundary are calculated as a function of grain size. The atomistic simulations show that the reduction of total elastic moduli with decreasing grain size is mainly due to a resulting larger grain-boundary atoms fraction, and that the total elastic moduli can be approximated by a simple weighted average of larger grains elastic moduli and a lower grain-boundary elastic moduli.


I. INTRODUCTION
Understanding the elastic properties of metals plays a key role in the analysis and modeling of the response to deformation under varying conditions.The local elastic properties in polycrystalline metals depend on the local structure, which varies from crystalline to grain boundaries and hetero-phase interfaces [1][2][3][4][5][6][7][8].The local variation of elastic properties plays a key role in tailoring of the macroscopic effective material properties, specifically in cases where grain refinement is used [6].Similarly local properties control strengthening in composite materials whether produced by additive manufacturing [9] or via self organized segregation (e.g. [10]) as well as in developing models for structure evolution [11].In order to develop an understanding of global response functions for composite materials, it is necessary to describe detailed distributions of local properties in addition to global averages.Since local elastic properties are not accessible experimentally for a wide range of systems, it is beneficial to study surrogate model systems, which allow reliable numerical evaluation of these.Such simulations play a key role in developing effective models for composite materials.
In Ref.
[12] we demonstrated the feasibility, robustness, and accuracy of calculating bulk elastic constants using molecular dynamics simulations in the NVT ensemble (where the number of particles, volume and temperature are kept constant).This method is generalized in this work to include the calculation of local elastic properties.The main advantage of this approach is that all components of the elasticity tensor are obtained in a single consistent moleculardynamics simulation, as opposed to the standard explicit deformation method, which requires several simulations under different deformations and is limited to evaluating global average values.Using NVT local evaluators, calculations are inherently local and immediate, so localized averaging leading to spatial and temporal distributions is trivial.In addition, while special care needs to be taken in the direct drive method to avoid strain rate and defect formation effects, the NVT local evaluation does not rely on forcing deformation.It thus is not affected by these potential problems.Moreover, in the standard deformation method, a non-negligible deformation of the simulation box is required in order to create a measurable variation in local stress.In inhomogeneous systems where elastic properties vary locally, such deformations lead to strain localization, which in many cases results in local yield [7, 13] that prevents the calculation of elastic moduli.
Previously, direct calculations of local elastic constants by atomistic simulations were performed in Refs.[14,15] for copper and gold, in Refs [16,17] for amorphous polymeric glasses, in Ref. [18] for a lipid bilayer and in Ref. [19] for ionic liquids.In Ref. [20] an experimental correlative approach was used to extract the local elastic properties of titanium.
In this paper, we present and analyze the distribution of local elastic constants in nanocrystalline copper and tantalum, as calculated using molecular dynamics simulations employing realistic many-body potentials.We extract the distributions of elastic moduli within polycrystalline systems and study their dependence on grain size in the range of 5-20 nm.In addition, we study the grain-size dependence of the average elastic moduli of the total system, in grains and grain boundary, and show explicitly that average global values agree with mean field models that were used in the literature [21][22][23][24][25][26].We note that existing models for the calculation of grain-boundary elastic properties of nanoscrystalline metals, usually rely on such effective mean field models, which assume a grain size independent grain-boundary moduli [27].Such methods only allow the calculation of global averages.The method which we employ in this work, gives a direct evaluation of local elastic properties, without relying on some of the major assumptions which were previously made -i.e. the description of the elastic properties of the grain boundary atoms using a single average value, or the independence of local properties on interface curvature.The method enables the calculation of new observables such as the distributions of elastic properties, which can play a crucial role in developing stochastic models addressing the response under nondeterministic drive conditions [11,28].

II. LOCAL ELASTICITY
In this section, we outline the method for the calculation of local elastic constants using molecular dynamics simulations in the NVT ensemble.
It was shown in previous works [12,15,16,18,[29][30][31][32][33][34], that in the NVT ensemble, the elasticity tensor can be written as a sum of three contributions: (i) a purely configurational part known as the Born term, which is given by a canonical average of the second order derivative of the potential energy with respect to the Lagrangian strain tensor, (ii) a stress fluctuation term that vanishes at zero temperature and (iii) a kinetic ideal gas contribution, which also vanishes at zero temperature.As a result, the low-temperature limit for the total elasticity tensor can be written on a per-atom basis [34,35], in the form: where α, β, γ, δ are the directional tensor indices, ⟨•⟩ represents ensemble average, C B i,αβγδ is the local, peratom Born elasticity tensor, V is the total system volume and V i is the volume of atom i, defined such that V = i V i .Following Refs.[34,35], we define the local volume of a specific atom as the volume of the Voronoi cell associated with it.In this work, calculations were performed for copper and tantalum, modeled by embedded-atom-model (EAM) potentials [36].These potentials are defined by a pair potential function v = v (r), an embedding function F = F (ρ) and a local density function ρ = ρ (r), so that the the potential energy takes the form: where r i is the position of atom i, r ij = r i − r j and ρ i = j̸ =i ρ (r ij ) is the density around atom i.It can be shown that for an EAM potential, the per-atom Born elasticity tensor takes the form [12,15,30,34,37,38]: where: and: Various local elastic moduli can be obtained directly from the local elasticity tensor (3), by employing local Voigt averages [39][40][41], which results in the following expressions for the local Bulk modulus: the local Shear modulus: the local Young's modulus: and the local Poisson's ratio: The Voigt notation for tensor indices was used in equations ( 6)-(7).

III. RESULTS
Simulations were performed using the state-of-the-art molecular dynamics code LAMMPS (Large-scale Atomic Molecular Massively Parallel Simulator) [42], to which we added a parallel implementation of the calculation of the local elasticity tensor for EAM potentials (equations (3)-( 5)).Results are reported for nanocrystalline copper and tantalum with grain sizes of d = 5, 8, 10, 12, 15, 18, 20nm and a corresponding number of atoms in the range of 2 × 10 5 − 1.5 × 10 7 .Copper was simulated using EAM potential by Mishin et al. [43] for an FCC lattice structure with a = 3.615A and tantalum was simulated using EAM potential by Ravelo et al. [44] for a BCC lattice with a = 3.304A.Initial configurations were generated by randomly assigning FCC and BCC crystallites on a grain-size super lattice and defining interfaces between grains using Voronoi tessellation.These initial configurations were relaxed using varying box-size NPT molecular dynamics simulations with a target pressure of zero for a duration of 200ps (with a 2fs time-step), as shown in Figure 1.It is evident that as the system reaches the state of equilibrium, the instantaneous pressure oscillates indefinitely with a well defined frequency and amplitude, while the cumulative temporal average pressure decays to zero.This behaviour of the pressure fluctuations is expected from a molecular dynamics simulation in the NVT ensemble (as shown in detail in Ref. [12], and in the references therein).
The final local pressure distributions in grain and grain boundaries are shown in Figure 5. Final relaxed configurations of several simulation boxes, visualized with Ovito [45] are shown in Figure 2, where atoms color are assigned according to the local average centrosymmetry parameter [46,47], distinctively separating grain and grain-boundary atoms [48].
In Figure 3, planar slices are shown with atoms colored according to the per-atom shear modulus (Eq.( 7)).
It is evident that the shear modulus varies considerably within grain-boundary regions in contrast with a significantly lower variation in the inter-grain regions.This behavior is observed for all elastic moduli (defined in equations ( 6)-( 9)), and is expected due to the contrast between the amorphous and crystalline atomic arrangements.
Histograms for the per-atom shear modulus of grain and grain-boundary atoms are shown in Figure 4, for various grain sizes.The distributions for intra-grain moduli have a width to average ratio of ∼ 5%, while the grain-boundary moduli distributions have a ratio close to unity.While this variation is expected, we further examined the dependence of these distributions on grain size.Here, we found that while intra-grain moduli distributions became sharper (lower width over mean) with increasing grain size, such variation was not observed in grain-boundary moduli distributions.Similar behavior is observed for the stress distributions as seen in Figure 5.We note that even though the average grainboundary shear modulus is lower than the average grain shear modulus (as expected), it is seen in Figure 4 that some grain-boundary local values are larger than local grain values.This is a result of the highly frustrated local amorphous configuration, which results in some grain boundary local values that are higher than those of the ordered phase.
We now discuss the uncertainty in the results.We expect the main variation to be due to the structure of interfaces and triple junctions formed during the Voronoi tessellation process.In order to quantify this uncertainty, we divided the simulation box into 8 different octants and compared the resulting moduli distributions in the various octants, which contain different grains, interfaces and junctions.We found that the different distributions are very similar, with the average value varying up to 2%, which can serve as a measure of the uncertainty of the average elasticity values (which are shown below in Fig. 8).This variation does not depend strongly on grain size.This is demonstrated in Fig. 6, where we compare the grain and grain boundary shear modulus distributions in 8 different octants and in the entire simulation box, for d = 10nm nanocrystlline copper.
In Figure 7, the atomic fraction of grain and grainboundary atoms are plotted as a function of grain size.It is evident that the grain-boundary atoms fraction has a strong decreasing dependence on grain size, as expected [48,49].In fact, the number of grain-boundary atoms must scale as the inverse of grain size; that is, if we denote by x g and x gb , the grain and grain-boundary atomic fractions, respectively, then: where d 0 is a constant length scale that depends on the lattice symmetry, and given by d 0 ≈ n gb ng ∆ gb where ∆ gb is the typical grain-boundary width and n gb , n g are, respectively, the grain and grain-boundary atom density.These results can be derived from a simplistic spherical grain-boundary model for which x gb = 4πd 2 ∆ gb n gb 4 3 πd 3 ng .From the data presented in Figure 7, we find that for copper d 0 ≈ 1.5nm and for tantalum d 0 ≈ 1.7nm.The scaling law (10) is also plotted in Figure 7, which shows that it perfectly matches the simulations' values.Finally, we also note that the decrease of grain-boundary atomic fraction with grain size d, can also be seen visually in Figures 2-3. Figure 8 presents the grain size dependence of various elastic moduli, namely, the shear modulus G, Young's modulus E and bulk modulus B, as well as the Poisson's ratio v.These moduli are plotted separately for grain and grain-boundary atoms and for the entire system (that is, all atoms in the simulations), resulting in the respective total bulk moduli.
It is evident that the total elastic moduli increases with grain size, in accordance with the results of Refs.[23,24,26,50].In order to analyze this behavior, we write the total elastic moduli as a sum of grain and grain-boundary moduli: where M represents one of the considered moduli G, E, B, v and M total , M g and M gb are, respectively, the total, grain and grain-boundary elastic moduli.The results in Figure 8 indicate that the grains have elastic moduli M g that are about 10% larger than the grain-boundary values M gb .This is to be expected since it is well-known that the grain boundary has an amorphous structure [1, 3, 4], a result that is consistent with our simulations and with Figure 4 which shows wide elastic modulus distributions for grain-boundary atoms.Moreover, it is evident that M g (d) and M gb (d) both have a relatively weak dependence on the grain size d.Therefore, we define a "mean-field" model, assuming grain-size independent grain and grain-boundary elastic moduli M g (d) ≈ M g , M gb (d) ≈ M gb , where M represents arbitrary average of the elastic modulus over grain size d.We take the arithmetic average over the grain sizes considered.Using the 1/d scaling of the grain-boundary atoms fraction (Eq.( 10)), we obtain a simple approximate mean-field model for the exact elastic moduli in Eq. ( 11), as a function of grain size: The results of this simple mean field model are also shown in Figure 8, offering a very good qualitative and even quantitative agreement with the exact results obtained from the simulations.Therefore, we deduce that the grain size dependence of the total elastic moduli M tot (d), is mainly due to the strong grainsize dependence of the grain and grain-boundary atoms fraction (via.Eq. ( 10)), while the effect of the grain size dependence of the grain and grain-boundary elastic moduli is only a second order effect.

IV. SUMMARY
Local elastic moduli of nanocrystalline copper and tantalum were calculated using molecular dynamics simulations.The resulting moduli distributions within the polycrystals were calculated and analyzed as a function of grain size in the range of 5-20 nm.The shear modulus distribution in the amorphous grain boundary is extended and grain-size independent, while the shear modulus distribution for atoms within the grains is significantly narrower and becomes sharper with increasing grain size.Even though local properties do not depend on grain size, global average elastic moduli do show grain size dependence in accordance with previous observations [16,17,51,52].This effect is shown explicitly to be denominated by the variation in relative numbers of crystalline versus grain-boundary atoms.We show for the simulated system that the average global elastic moduli can be approximated by a weighted average of the larger grains elastic moduli and the lower grain-boundary elastic moduli, an approximation which was assumed in several previous works [21,23,50,53,54].Moreover, by employing the simple inverse relation between the grain-boundary atoms fraction and grain size d, we suggest a simple but accurate mean-field formula for the grain size dependence of various bulk atomic moduli (shear modulus, Young's modulus, bulk modulus and Poisson's ratio).* menahem.krief@mail.huji.ac.il [1] Herbert Gleiter.Nanostructured materials: basic concepts and microstructure.Acta materialia, 48(1):1-29, 2000.

FIG. 1 .
FIG.1.Demonstration of the equilibration process for nanocrystalline tantalum with grain size d = 100nm.The box size is relaxed in the NPT ensemble, until the total pressure is zero.Shown are the total instantaneous pressure (blue solid line, left y-axis) and the average pressure (red dashed line, right y-axis), which is calculated over a window of 10 4 times steps (with a time step of 2fs).

FIG. 2 .
FIG. 2. Atomistic visualizations of simulated nanocrystalline copper (FCC, upper panes) and tantalum (BCC, lower panes), with grain sizes d = 8nm (left panes), 15nm (middle panes) and 20nm (right panes).The corresponding box size L is also listed in each pane.Color indicates the local average centrosymmetry parameter.It is evident that the number of grains is kept constant for different grains sizes.

FIG. 3 .
FIG.3.Atomistic view of (111) planar slices of simulated nanocrystalline copper (upper panes) and tantalum (lower panes), with grain sizes d = 8nm (left panes), 15nm (middle panes) and 20nm (right panes).The corresponding simulation box size L is also listed in each pane.The spatial distribution of the local shear modulus G, is indicated by color.The lack of spatial correlation between the widely distributed values within the grain boundary atoms is evident.In contrast, deviations from the average in the grain are significantly smaller, as expected.It is also seen that the grain atoms have a larger value on average.These results are consistent with Fig.4.We note that the parallel lines in some of the grains are due to imgae shading, and are artificial.

FIG. 4 .
FIG. 4.Shear modulus probability densities for grain and grain-boundary atoms, in nanocrystalline copper (upper pane) and tantalum (lower pane) and for various grain sizes d (as listed in the legend).
FIG. 7.The fractional number of grain (in blue) and grain-boundary (in red) atoms as a function of grain size in nanocrystalline copper (upper pane) and tantalum (lower pane).The points are the atomic fractions in the simulations (see Figure2), and the dashed lines are obtained from a 1/d fit (according to the scaling relation (10), as discussed in the text).

FIG. 8 .
FIG.8.Elastic moduli for grain-boundary (red circles), grain (blue squares) and all (black diamonds) atoms, as a function of grain size, for nanocrystalline copper (upper figures) and tantalum (lower figures).From left to right, the shear modulus G, Young's modulus E, bulk modulus B, and Poisson's ratio v are shown.The dashed black line represents the mean-field model with 1/d scaling of the grain boundary contribution (Eq.(12)) with d0 ≈ 1.5nm for copper and d0 ≈ 1.7nm for tantalum).