Incipience of Plastic Flow in Aluminum with Nanopores: Molecular Dynamics and Machine-Learning-Based Description

Abstract

Incipience of plastic flow in nanoporous metals under tension is an important point for the development of mechanical models of dynamic (spall) fracture. Here we study axisymmetric deformation with tension of nanoporous aluminum with different shapes and sizes of nanopores by means of molecular dynamics (MD) simulations. Random deformation paths explore a sector of tensile loading in the deformation space. The obtained MD data are used to train an artificial neural network (ANN), which approximates both an elastic stress–strain relationship in the form of tensor equation of state and a nucleation strain distance function. This ANN allows us to describe the elastic stage of deformation and the transition to the plastic flow, while the following plastic deformation and growth of pores are described by means of a kinetic model of plasticity and fracture. The parameters of this plasticity and fracture model are identified by the statistical Bayesian approach, using MD curves as the training data set. The present research uses a machine-learning-based approximation of MD data to propose a possible framework for construction of mechanical models of spall fracture in metals.

1. Introduction

Dislocations are known to be the main defects of crystal structures responsible for the plasticity of metals [1,2]. If a crystal’s interior is initially devoid of dislocations, they have to be nucleated homogeneously inside perfect areas of crystal [3,4,5,6,7,8,9,10,11] or be emitted from other defects such as grain boundaries [12,13,14,15,16,17,18,19], phase inclusions [20,21,22,23], interfaces [24,25,26], and pores [27,28,29,30] in order to start the plastic flow. Even if a crystal is not absolutely devoid of dislocations, the pre-existing dislocation density can be not enough to accommodate either ultra high rate dynamic loading with the strain rate up to 10⁹ s⁻¹ [31,32,33,34,35] or small-scale localized loading at nanoindentation [36,37,38,39,40]. In the first case, the rate of plastic slip and the rate of multiplication of the pre-existing dislocations can be insufficient, and the growing shear stress can activate either homogeneous nucleation of dislocations or heterogeneous nucleation and emission. In the second case, local loading touches dislocation-free areas of crystal with a high probability. In both cases, the nucleation of dislocations triggers and controls the plastic flow, which motivates the study of this phenomenon.

Homogeneous nucleation of dislocations in perfect single crystals starts at relatively high elastic deformation: the elastic shear strain reaches 5–8% for aluminum [10] and copper [41,42]. This elastic shear makes the crystal lattice unstable to the formation of nano-scale dislocation loops, which starts the plastic flow. The instability is preceded by a nonlinear increase or even local decrease in shear stress with shear strain approaching the nucleation threshold. This nonlinear behavior leads to a decrease in effective shear modulus and, consequently, to a decrease in the barrier of dislocation nucleation [5,10]. Previous studies show substantial dependence of the dislocation nucleation threshold on the loading conditions. Density functional theory (DFT) calculations [43] show the influence of tri-axial and uniaxial stresses on the ultimate shear strength of copper and aluminum. A tension–compression asymmetry of dislocation nucleation was first reported in [3,4] and explained in [7]; classical molecular dynamics (MD) simulations were utilized in these papers for investigation of the dislocation nucleation. Classical MD allows one to consider multi-atom systems and examine a lot of configurations with relatively moderate computational cost in comparison with DFT. On the other hand, the MD simulations reveal a lot of atomistic details in the deformation process and are highly suitable for the problem of dislocation nucleation. Insufficiency of two stress components, the resolved shear stress (RSS) and the resolved normal stress (RNS), for determination of the nucleation condition is concluded in [44] based on the analyses of MD simulations. A 4-fold increase in the RSS (from about 1.5 GPa to about 6 GPa) is reported in [41] for copper as the pressure rises from the negative level of −10 GPa to the positive level of +50 GPa. Previously, we showed [41,42,45] that study of various deformation paths by means of MD simulations and generalization of the obtained results in the form of a trained artificial neural network (ANN) is an efficient approach to construct the dislocation nucleation criterion. The MD provides big data for training, while the ANN is an optimal tool for reproducing a complex functional dependence. Although we discuss the nucleation threshold as a certain condition for the plasticity incipience, nucleation itself is a stochastic process and depends on the strain rate. Concurrently, the strain rate dependence is weak and can be taken into account as proposed in [41] if we know the nucleation threshold at any fixed strain rate.

Various surface and volumetric defects substantially reduce the threshold of dislocation nucleation. This is due to the fact that these defects can act as stress concentrators; additionally, the atoms in their vicinity already have increased energy, which reduces the energy barrier of dislocation nucleation. Plasticity incipience around a pore in metal starts compacting this pore under the compressive loading of metal [30,46,47] and, vice versa, pores grow under tensile loading of metal [27,29]. The first process is essential for dynamic compaction of porous metals, while the second process is a stage of dynamic tensile (spall) fracture of metals. MD investigation of the influence of the shock wave direction on the dislocation nucleation and compaction of pores in nanoporous aluminum single crystals is reported in [47]. Compaction of nanoporous aluminum and the influence of nano-sized pores on the threshold of dislocation nucleation is investigated in [30] by means of MD simulations in combination with ANN training. Only compressive loading is considered in this paper, and only uniaxial and uniform tri-axial deformations among the possible deformation paths are analyzed. The study of the tensile loading of nanoporous metals along various deformation paths is insufficiently covered in the literature. This study is necessary for both a deep understanding of the dislocation nucleation and for further development of the mechanical models of dynamic (spall) fracture [48,49,50]. Activation of nanopore growth on pre-existing volumetric defects instead of homogeneous nucleation of nanopores is the main mechanism of fracture initiation at moderate strain rates as shown in [49,50]. Therefore, the influence of nanopores on the plasticity incipience is even more relevant.

In this paper, we adopt the method of random deformation paths proposed in [51] in order to explore a sector of tensile loading in the deformation space. By means of MD simulations, we consider axisymmetric deformation of nanoporous aluminum with different shapes and sizes of nanopores. The obtained MD data are used to train an ANN approximating a tensor equation of state [42,45] of porous aluminum and a nucleation strain distance function [41,42]. This ANN allows us to describe the elastic stage of deformation and the transition to the plastic flow, while the following plastic deformation and growth of pores are described by means of a kinetic model of plasticity and fracture. The parameters of this fracture model are identified by means of a statistical Bayesian approach [50,51,52,53], using MD curves as the training data. The present research uses a machine learning-based approximation of MD data to propose a possible framework for construction of mechanical models of spall fracture in metals.

2. Materials and Methods

2.1. Molecular Dynamics Simulations

Molecular dynamics (MD) simulations are used to explore the plasticity incipience, as well as the elastic and plastic deformation behavior of nanoporous aluminum. The LAMMPS software package [54] (Sandia National Laboratories, Albuquerque, NM, USA) is used to perform MD simulations together with the widely used interatomic potential [55] for Al atoms. This potential belongs to the EAM type [56] and takes into account the interaction of atoms with a density of free electrons specific for metals. The used potential was previously verified in [10] for the problem of dislocation nucleation in a perfect single crystal of aluminum in comparison with the results of the first-principle calculations taken from [43,57]. In Section 2.3, we also compare MD results with the density functional theory (DFT) calculations for the hydrostatic tension of aluminum.

The basic MD system contains 500,000 atoms of aluminum and has dimensions of 50 × 50 × 50 lattice constants (about 20 × 20 × 20 nm³) for the FCC lattice of aluminum. We consider a perfect single crystal with the lattice directions [100], [010], and [001] coinciding with the $x$, $y$, and $z$ coordinate axes, respectively. Periodic boundary conditions applied to all boundaries of the MD system make it equivalent to a small representative volume element (RVE) of an infinite bulk material. Besides pore-free solid aluminum, we consider nanoporous samples with spherical pores, cylindrical pores (straight circular cylinders with diameter equal to the base height), and cubic pores. For each pore shape, three different equivalent sizes are analyzed: 4, 8, and 16 nm. The equivalent size is the diameter of a spherical pore of the same volume. The considered equivalent sizes provide the initial porosities of about 0.42%, 3.3%, and 27%, respectively. Thus, a wide range of porosities is considered. The atoms inside the pore are deleted from the MD system at its construction by means of the “delete_atoms” command of LAMMPS, such that the lattice around the pore remains perfect and the number of atoms decreases proportionally to the initial porosity.

Before deformation, the MD system is thermalized at 10 ps at the test temperature and zero stresses in an NPT ensemble by controlling of the number of atoms, stresses, and temperature. Both room temperature (300 K) and an elevated temperature of 700 K are considered. Thereafter, the MD system is axisymmetrically deformed in an NVT ensemble by means of the “deform” command in LAMMPS, which rescales the coordinates of all atoms and avoids mechanical waves and parasitic vibrations in the system. This deformation is physically equivalent to the expansion due to inertia.

An axisymmetric deformation is determined by two independent engineering strains, ${\epsilon }_1$ and ${\epsilon }_2$, as follows:

$${\epsilon }_1=\frac{L_x-L_{x0}}{L_{x0}}=F_{11}-1,{\epsilon }_2=\frac{L_y-L_{y0}}{L_{y0}}=\frac{L_z-L_{z0}}{L_{z0}}=F_{22}-1,$$

(1)

where $L_x$, $L_y$, and $L_z$ are the current lengths of the RVE along the corresponding axes; $L_{x0}$, $L_{y0}$, and $L_{z0}$ are the initial lengths; $F_{11}$ and $F_{22}=F_{33}$ are the diagonal components of the total deformation gradient, whereas non-diagonal components are zero. The engineering strain rates, ${\dot{\epsilon }}_1$ and ${\dot{\epsilon }}_2$, are chosen randomly, and they are constrained as in [42] the following requirement:

$$\left|{\dot{\epsilon }}_1\right|+2\left|{\dot{\epsilon }}_2\right|=\dot{\epsilon }={10}^9{\mathrm{s}}^{-1},$$

(2)

where $\dot{\epsilon }$ is the total engineering strain rate. We characterize a deformed state using the Green–Lagrange strain tensor, which has two independent diagonal components:

$$E_{11}=\frac{1}{2}\left[F_{11}^2-1\right]={\epsilon }_1+\frac{1}{2}{\epsilon }_1^2, E_{22}=E_{33}=\frac{1}{2}\left[F_{22}^2-1\right]={\epsilon }_2+\frac{1}{2}{\epsilon }_2^2$$

(3)

and zero non-diagonal components in the considered case. The deformation paths in the $\left\{E_{11},E_{22}\right\}$ plane are shown in Figure 1. We study 20 random paths for nanoporous aluminum and 30 paths for pore-free aluminum. The total number of MD simulations reaches 420, including 3 shapes and 3 sizes of pores, 2 temperatures and 20 deformation paths for nanoporous samples, and 2 temperatures with 30 deformation paths for pore-free ones.

The Nose–Hoover thermostat and barostat [58] are used to control the temperature and stresses when necessary. The pressure averaged over the metal volume including pores is calculated by means of the Virial theorem [59]. The OVITO program [60] is used for both visualization and analysis of the atomic configurations. The dislocation extraction algorithm (DXA) [61] is used to detect dislocations in the crystal lattice, which allows us to find out the dislocation nucleation event and to trace further evolution of the dislocation density. The construct surface mesh algorithm [62] is applied to calculate the current void volume and porosity. For visualization, we use ranking of atoms in accordance with the centrosymmetry parameter [63], which shows the degree of deviation of the atom’s environment from an ideal crystal lattice.

Figure 2, Figure 3 and Figure 4 show several examples of considered MD systems and the processes of dislocation nucleation and emission from the pore surface. The plasticity incipience begins with the formation of half-loops of the leading partial Shockley dislocations, leaving the stacking faults; however, split perfect dislocations then predominate. One can see that the deformation path and the pore shape influence the dislocation nucleation.

2.2. Processing of MD Data and Approximation by Artificial Neural Network

Processing of MD data gives us the evolution of the volume-averaged longitudinal stress ${\sigma }_{11}$, the volume-averaged transverse stress ${\sigma }_{22}={\sigma }_{33}$, the specific energy $U$, the dislocation density ${\rho }_{\mathrm{D}}$, and porosity $\phi$ in the RVE. These data are used to construct a tensor equation of state in the elastic range of deformations and an indicator of plasticity incipience in the form of a trained artificial neural network (ANN). In some cases, dislocations can nucleate on the pore surface and then disappear or oscillate for some time before the start of the obvious plastic stage. In order to disregard these preliminary processes, we treat the formation of the total dislocation length above a conditional level of 3 nm in the RVE as the plasticity incipience. Let us denote as ${\epsilon }_{\mathrm{n}}$ the engineering strain at the moment of plasticity incipience. It is proposed in [41] to use the difference $\epsilon -{\epsilon }_{\mathrm{n}}$ between the current strain $\epsilon$ along a certain deformation path and the nucleation strain ${\epsilon }_{\mathrm{n}}$ as the nucleation strain distance function $Q$. This function becomes positive beyond the threshold of dislocation nucleation at the considered strain rate, and its slope allows one to take into account the strain rate dependence of nucleation as shown in [41]. The definition $Q=\epsilon -{\epsilon }_{\mathrm{n}}$ leads to an ambiguity in the training data at $\epsilon =0$, because it gives different values $Q=-{\epsilon }_{\mathrm{n}}$ for different deformation paths. In order to eliminate this ambiguity, we redefine the nucleation strain distance function as follows:

$$Q=\max \left\{\epsilon -{\epsilon }_{\mathrm{n}}, -0.05\right\}$$

(4)

where an artificial minimum of −0.05 provides uniqueness of $Q$ at small deformations, whereas it does not influence the behavior of $Q$ near the nucleation threshold. Only this behavior is required for application of the nucleation strain distance function.

Plasticity influences the evolution of stresses and energy beyond the nucleation threshold. The stresses become dependent on the previous plastic deformation and cease to be a single-valued function of the strain tensor. We intend to construct a tensor equation of state as the dependence of stresses and energy on the elastic strains. Accordingly, we distinguish the elastic parts of the deformation curves and approximate the dependencies of stresses and energy along each deformation path as third-order polynomials of the engineering strain. The coefficients of polynomials are calculated by the least square method. Thereafter, these polynomials are used to prepare the training data for the ANN in both the elastic domain and beyond the nucleation threshold up to $1.5{\epsilon }_{\mathrm{n}}$. We therefore extrapolate the elastic behavior by means of these polynomials at $\epsilon >{\epsilon }_{\mathrm{n}}$. Using polynomial approximation for each separate deformation curve is simple. An attempt to apply a similar approximation of the results of all 420 MD runs simultaneously is doomed to fail, whereas the application of the ANN is an efficient method.

The axial stress ${\sigma }_{11}$, the longitudinal stress ${\sigma }_{22}$, the specific energy $U$, and the nucleation strain distance function $Q$ form the output vector $Y=\left\{{\sigma }_{11},{\sigma }_{22},U,Q\right\}$ of the ANN. The form factor $\kappa$, the initial porosity ${\phi }_0$, the longitudinal $E_{11}$ and transverse $E_{22}$ components of the Green–Lagrange strain tensor (see Equation (3)), and the temperature $T$ form the input vector $X=\left\{\kappa ,{\phi }_0,E_{11},E_{22},T\right\}$ of the ANN. The form factor $\kappa$ is defined as the proportionality coefficient between the volume of a pore $V_{\mathrm{p}}$ and the cube of its linear dimension $a$ (diameter, edge length, or height of the cylinder): $V_{\mathrm{p}}=\kappa a^3$. It has the value of $\kappa =\pi /6$ for spherical pores, $\kappa =\pi /4$ for cylindrical pores, and $\kappa =1$ for cubic pores. The introduced above equivalent size $d$ of a pore can be calculated as $d=a\sqrt[3]{6\kappa /\pi }$. The pairs of input and output vectors $\left\{X^{\beta },Y^{\beta }; \beta \in \left[1,\mathsf{{\rm B}}\right]\right\}$ following from the processing of MD simulations form the training data set. The training data are defined in the range $\epsilon \in \left[0, 1.5{\epsilon }_{\mathrm{n}}\right]$ with a step of $0.02{\epsilon }_{\mathrm{n}}$ along each deformation path. The data for pore-free aluminum are repeated three times with zero initial porosity and all three values of the form factor because systems with all shapes of pores must tend to the case of pore-free aluminum at zero porosity. The total number $\mathsf{{\rm B}}$ of the input–output pairs in the prepared training data set reaches about 40,000. The training data are collected in a worksheet called “Training_Data” in the file “NPA22.xls” (Supplementary Materials).

The feedforward ANN we use is a complex superposition of simple transfer functions. The structure of the ANN consists of an input layer containing 5 neurons according to the dimension of the input vector $X$, 6 hidden layers with 20 neurons in each layer, and an output layer containing 4 neurons according to the dimension of the output vector $Y$. The “Leaky ReLU” transfer function is used for all neurons of hidden layers, whereas the “Sigmoid” transfer function is used for the output layer. The ANN maps $X\to Y$, and the training consists of a minimization of the discrepancy of the ANN results on the training data by means of a successive random variation of the weights and biases of neurons. The minimization procedure can be expressed as a set of matrix operations commonly known as the back propagation algorithm [64,65]. The procedure that we use for ANN training is described in detail elsewhere [10,41,45]. Application of this training procedure, realized by in-house single-thread FORTRAN computer code, allows us to reach an average error of 0.35% and a maximum error of 6.5% in 5 h of training. The ANN is ultimately designed to model continuum mechanics within our own software packages, which encourages use of our own code to train the ANN to achieve full compatibility. The obtained weights and biases of neurons are collected as a worksheet named “Al.TEOS1.ANNp” in the file “NPA22.xls” (Supplementary Materials); the structure of this data file is described in [10,41,42]. Figure 5 illustrates the quality of ANN training using an example of longitudinal stress ${\sigma }_{11}$. The ANN results in comparison with the training data from the MD are collected in the worksheet named “ANN_Results” in the file “NPA22.xls” (Supplementary Materials).

2.3. Density Functional Theory Calculations

Precision of MD simulations greatly depends on the quality of interatomic potential. Calculations based on the density functional theory (DFT) are more reliable because the DFT approximately solves the quantum mechanical problem for atomic electrons [66]. The DFT employs a number of approximations, but these calculations are usually considered as first principles or ab initio ones. DFT calculations of large atomic systems, such as the RVE of porous metal, and accounting of finite temperatures are still challenging problems, but the study of various deformed states of pore-free metal at 0 K is a promising application of the DFT in exploring the tensor cold curve. Such DFT-based cold curves can be used for construction of a tensor equation of state in future works. We use DFT calculations in order to verify the trained ANN for the tensor equation of state of nanoporous aluminum and, consequently, the interatomic potential. We calculate a number of deformed states with hydrostatic ($E_{11}=E_{22}$) tension and compare the obtained pressure-strain dependence with the results of the ANN in Figure 6.

Quantum ESPRESSO [67,68] software, the GPU-optimized version [69], is applied for DFT calculations. A one-atom system using FCC symmetry is considered. The calculations are carried out using the self-consistent field (SCF) method on plane waves (PW) with the Perdew–Burke–Ernzerhof (PBE) exchange–correlation energy functional [70] until convergence of the total energy expressed in Ry at eight decimal places. The influence of ions on the conduction electrons is taken into account by means of the Dal Corso pseudopotential for Al from pslibrary [71] with Methfessel–Paxton smearing [72]. The K-space grid is 20 × 20 × 20 points. The cut-off energy for wave functions is 29 Ry, and the cut-off energy for electron density is 290 Ry.

A close fit between the DFT-calculated cold curve and the ANN for 300 K is observed for strains below $0.8 {\epsilon }_{\mathrm{n}}$ which verifies both the trained ANN and the used interatomic potential [55]. The following divergence can be explained by the influence of thermal fluctuations leading to the nucleation of dislocations and voids in multi-atom systems at room temperature. Note that the ANN reproduces the polynomial extrapolation of the pressure curve beyond the nucleation threshold. Also note that there is no pure hydrostatic loading among the random MD deformation paths used for ANN training.

2.4. Plasticity and Fracture Model

Description of plasticity is based on the modified Maxwell relaxation model. The used modification is proposed in [73] for metals and takes into account a threshold of plastic flow determined by the yield strength, as well as dependence of the yield strength and relaxation time on the dislocation density. This dislocation density dynamically changes in accordance with the kinetics equations that complement the model. The model of porosity variation is formulated in [30] for various shapes of pores and for arbitrary tri-axial deformation; this model develops the previous approach [27], which distinguishes the current and plastically stabilized dimensions of pores, as well as a gradual dislocation-driven plastic change of the stabilized dimensions towards the current ones. The combined model of porosity variation and macroscopic plasticity is formulated in [30] in terms of finite deformations and applied to the problem of dynamic compaction of nanoporous aluminum. Here we apply the model to the opposite case of fracture with increasing porosity over time. In addition, we replace the linear Hooke’s law used in [30] with the nonlinear tensor equation of state in the form of the ANN and use the ANN-based description of the dislocation nucleation as in other papers [42] devoted to the plastic deformation of solid copper. Therefore, the used model is a development of the previous ones, and we formulate it here in detail, which is necessary for reproducing the model. We restrict ourselves to the axisymmetric case of deformation.

The total (macroscopic) deformation gradient is determined as derivatives of the current coordinates over initial ones and has two independent diagonal components, $F_{11}$ and $F_{22}=F_{33}$, whereas all non-diagonal components are zero for axisymmetric deformation. In the considered case, $F_{11}$ and $F_{22}$ are completely determined by the dynamics of tension of the RVE (see Equation (1)). A multiplicative decomposition [74,75,76,77] is used, which supposes intermediate mappings of coordinates responsible for different physical processes, such as plasticity characterized by $F_{11,22}^{\mathrm{p}}$ and damage (porosity) development characterized by $F_{11,22}^{\mathrm{d}}$. The elastic part of the deformation gradient and the corresponding Green–Lagrange strain tensor components $E_{11,22}$ can be calculated as follows [30]:

$$F_{11}^{\mathrm{e}}={\left(F_{11}^{\mathrm{p}}{F}_{11}^{\mathrm{d}}\right)}^{-1}{F}_{11}, F_{22}^{\mathrm{e}}={\left(F_{22}^{\mathrm{p}}{F}_{22}^{\mathrm{d}}\right)}^{-1}{F}_{22},$$

(5)

$$E_{11}=\frac{1}{2}[{\left(F_{11}^{\mathrm{e}}\right)}^2-1], E_{22}=\frac{1}{2}[{\left(F_{22}^{\mathrm{e}}\right)}^2-1].$$

(6)

The volume-averaged stresses ${\sigma }_{11,22}$ are calculated by means of the trained ANN with the strain tensor components as inputs as discussed above in Section 2.2. The solid phase occupies the part $\left(1-\phi \right)$ of the volume of porous material; therefore, the stresses in the solid phase can be estimated as follows:

$${\sigma }_{11}^{\mathrm{s}}=\frac{{\sigma }_{11}}{1-\phi }, {\sigma }_{22}^{\mathrm{s}}=\frac{{\sigma }_{22}}{1-\phi },$$

(7)

where ${\sigma }_{11}$ and ${\sigma }_{22}$ are the volume-averaged stresses, which are the outputs of the ANN.

Considering a pore characterized by a longitudinal size $a_1$ and transverse size $a_2$ with a volume $V_{\mathrm{p}}=\kappa a_1{a}_2^2$, where $\kappa$ is the form factor introduced above, the following kinematics equations for the damage deformation gradient are derived in [30]:

$$\begin{gathered} {\dot{F}}_{11}^{\mathrm{d}}=\frac{{\dot{F}}_{11}}{F_{11}^{\mathrm{e}}{F}_{11}^{\mathrm{p}}}\left[1-\frac{F_{11}^{\mathrm{d}}{\left(F_{22}^{\mathrm{d}}\right)}^2}{\left(1-{\phi }_0\right)}\right]+\frac{\kappa {\dot{a}}_1{a}_2^2}{F_{11}^{\mathrm{e}}{F}_{11}^{\mathrm{p}}{\left(F_{22}^{\mathrm{e}}{F}_{22}^{\mathrm{p}}\right)}^2}\frac{F_{11}^{\mathrm{d}}}{V_0\left(1-{\phi }_0\right)}, \\ {\dot{F}}_{22}^{\mathrm{d}}=\frac{{\dot{F}}_{22}}{F_{22}^{\mathrm{e}}{F}_{22}^{\mathrm{p}}}\left[1-\frac{F_{11}^{\mathrm{d}}{\left(F_{22}^{\mathrm{d}}\right)}^2}{\left(1-{\phi }_0\right)}\right]+\frac{\kappa a_1{\dot{a}}_2{a}_2}{F_{11}^{\mathrm{e}}{F}_{11}^{\mathrm{p}}{\left(F_{22}^{\mathrm{e}}{F}_{22}^{\mathrm{p}}\right)}^2}\frac{F_{22}^{\mathrm{d}}}{V_0\left(1-{\phi }_0\right)}, \end{gathered}$$

(8)

where ${\phi }_0$ is the initial porosity and $V_0$ is the initial volume of the RVE. The initial values are $F_{11}^{\mathrm{d}}=F_{22}^{\mathrm{d}}=1$.

The plastic deformation gradient can be written as follows:

$$F_{11}^{\mathrm{p}}=\exp \left(w\right), F_{22}^{\mathrm{p}}=\exp \left(-\frac{1}{2}w\right),$$

(9)

where $w$ is the accumulated plastic strain as the integral of the plastic strain rate $\dot{w}$ over time. The plastic strain rate $\dot{w}$ is calculated in accordance with the modified Maxwell relaxation model [73]:

$$\dot{w}={\eta }^{-1}\left({\tau }^{\mathrm{s}}-\frac{1}{2}Y·\mathrm{sign}\left({\tau }^{\mathrm{s}}\right)\right)·\mathsf{\Theta }\left(\left|{\tau }^{\mathrm{s}}\right|-\frac{1}{2}Y\right),$$

(10)

where ${\tau }^{\mathrm{s}}=\left({\sigma }_{11}^{\mathrm{s}}-{\sigma }_{22}^{\mathrm{s}}\right)/2$ is the average shear stress in the solid phase, which causes the average plasticity; $\mathsf{\Theta }(·)$ is the Heaviside step function. The static yield strength $Y$ is subject to strain hardening expressed in the form of the Taylor hardening law:

$$Y=\alpha Gb\sqrt{{\rho }_{\mathrm{D}}},$$

(11)

where $G$ is the shear modulus of the solid phase of aluminum, $b$ is the modulus of the Burgers vector of dislocations, $\alpha$ is the hardening coefficient, and ${\rho }_{\mathrm{D}}$ is the dislocation density. The viscosity coefficient $\eta$ is proportional to the relaxation time and inversely proportional to ${\rho }_{\mathrm{D}}$ as it is shown in [78]:

$${\eta }^{-1}=\frac{b{\rho }_{\mathrm{D}}^2}{4B},$$

(12)

where $B$ is the dislocation friction coefficient, which can be taken from the MD simulations [79] for aluminum.

Equation (8) must be complemented by the equations of variation of the pore sizes. Similar to [27,30], we distinguish the current sizes $a_{1,2}$ and the plastically stabilized ones $a_{1,2}^{\mathrm{p}}$, which are connected by means of the following relationship:

$$a_{1,2}=\left(a_{1,2}^{\mathrm{p}}-\frac{\gamma }{G}\right){\left[1+P\left(\frac{1}{4G}+\frac{1}{3K}\right)\right]}^{-1}$$

(13)

if one neglects the elastic vibrations of $a_{1,2}$. In Equation (13), $\gamma$ is the surface tension coefficient, $K$ is the bulk elastic modulus of the solid phase of aluminum, and $P=-\left({\sigma }_{11}+2{\sigma }_{22}\right)/3$ is the volume-averaged pressure. The difference between the current and plastically stabilized dimensions determines the elastic deformation around the pore and the corresponding shear stresses ${\tau }_{1,2}$ [27,30]:

$${\tau }_{1,2}=3G\left(1+\frac{P}{3K}-\frac{G{a}_{1,2}^{\mathrm{p}}}{G{a}_{1,2}^{\mathrm{p}}-\gamma }\left[1+P\left(\frac{1}{4G}+\frac{1}{3K}\right)\right]\right).$$

(14)

In the pore vicinity, the local shear stresses ${\tau }_{1,2}$ cause the plastic flows, with the plastic strain rates calculated similar to Equation (10):

$${\dot{w}}_{1,2}={\eta }^{-1}\left({\tau }_{1,2}-\frac{1}{2}Y·\mathrm{sign}\left({\tau }_{1,2}\right)\right)·\mathsf{\Theta }\left(\left|{\tau }_{1,2}\right|-\frac{1}{2}Y\right).$$

(15)

These plastic flows change the corresponding plastically stabilized dimensions $a_{1,2}^{\mathrm{p}}$ as follows:

$${\dot{a}}_{1,2}^{\mathrm{p}}=b{\dot{w}}_{1,2}.$$

(16)

The last portion of equations describes the dislocation kinetics. Variation of the dislocation density is described by a balance equation:

$${\dot{\rho }}_{\mathrm{D}}={\mathsf{{\rm P}}}_{\mathrm{n}}+{\mathsf{{\rm P}}}_{\mathrm{m}}-{\mathsf{{\rm P}}}_{\mathrm{a}},$$

(17)

where ${\mathsf{{\rm P}}}_{\mathrm{n}}$ is the nucleation rate, ${\mathsf{{\rm P}}}_{\mathrm{m}}$ is the multiplication rate, and ${\mathsf{{\rm P}}}_{\mathrm{a}}$ is the annihilation rate. The multiplication and annihilation rates are taken in the form proposed in [30]:

$${\mathsf{{\rm P}}}_{\mathrm{m}}=\frac{1}{{\epsilon }_{\mathrm{m}}}\left[{\tau }_{\mathrm{s}}\dot{w}+\frac{\phi }{2}\left({\tau }_1{\dot{w}}_1+{\tau }_2{\dot{w}}_2\right)\right],$$

(18)

$${\mathsf{{\rm P}}}_{\mathrm{a}}=k_{\mathrm{a}}{\rho }_{\mathrm{D}}\left[\left|\dot{w}\right|+\frac{\phi }{2}\left(\left|{\dot{w}}_1\right|+\left|{\dot{w}}_2\right|\right)\right],$$

(19)

where ${\epsilon }_{\mathrm{m}}$ is the dislocation multiplication energy, which is an efficient energy of formation of a dislocation line of unit length, and $k_{\mathrm{a}}$ is the annihilation coefficient.

Similar to [42], the nucleation term ${\mathsf{{\rm P}}}_{\mathrm{n}}$ is modified by using the nucleation strain distance function $Q$ represented in the form of the trained ANN as discussed above in Section 2.2. The function $Q$ describes the threshold of dislocation nucleation, and the nucleation rate can be written as follows [42]:

$${\mathsf{{\rm P}}}_{\mathrm{n}}=\frac{2\pi c_{\mathrm{t}}\rho }{m_{\mathrm{Al}}}·\exp \left(-\frac{k_{\mathrm{n}}G{b}^2}{k_{\mathrm{B}}T}\right)·\mathsf{\Theta }\left(Q\right),$$

(20)

where $c_{\mathrm{t}}=\sqrt{G/\rho }$ is the transverse sound speed, $\rho$ is the mass density of the solid phase, $m_{\mathrm{Al}}$ is the mass of one aluminum atom, $k_{\mathrm{n}}$ is a dimensionless coefficient of the nucleation energy barrier $k_{\mathrm{n}}G{b}^2$, and $k_{\mathrm{B}}$ is Boltzmann’s constant.

The parameters of the model, the surface tension coefficient $\gamma$, the dislocation multiplication energy ${\epsilon }_{\mathrm{m}}$, the hardening coefficient $\alpha$, the coefficient of dislocation nucleation barrier $k_{\mathrm{n}}$, and the annihilation coefficient $k_{\mathrm{a}}$ are determined by means of the Bayesian algorithm with MD curves as the training data. Other required parameters entering Equations (11)–(20) are either taken from the ANN or determined directly for the used interatomic potential by separate MD simulations.

The Bayesian algorithm involves enumeration of a number of random sets of model parameters $\mathsf{A}=\left\{ \gamma ,{\epsilon }_{\mathrm{m}},\alpha ,k_{\mathrm{n}},k_{\mathrm{a}}\right\}$ with estimation of probability $\mathsf{\Pi }\left(\mathsf{A}\right)$ as a measure of correspondence of the model predictions $\left\{{\sigma }_{11},{\sigma }_{22},\phi ,{\rho }_{\mathrm{D}}\right\}$ and the training data $\left\{{\sigma }_{11}^{\mathrm{MD}},{\sigma }_{22}^{\mathrm{MD}},{\phi }^{\mathrm{MD}},{\rho }_{\mathrm{D}}^{\mathrm{MD}}\right\}$ from the MD:

$$\mathsf{\Pi }\left(\mathsf{A}\right)={\displaystyle \prod }_{n=1}^N{\displaystyle \prod }_{m=0}^M\exp \left(-0.001\left[{\left(\frac{{\sigma }_{11}-{\sigma }_{11}^{\mathrm{MD}}}{\mathsf{\Delta }{\sigma }_{11}^{\mathrm{MD}}}\right)}^2+{\left(\frac{{\sigma }_{22}-{\sigma }_{22}^{\mathrm{MD}}}{\mathsf{\Delta }{\sigma }_{22}^{\mathrm{MD}}}\right)}^2+{\left(\frac{\phi -{\phi }^{\mathrm{MD}}}{\mathsf{\Delta }{\phi }^{\mathrm{MD}}}\right)}^2+{\left(\frac{\ln \left(1+{\rho }_{\mathrm{D}}\right)-\ln \left(1+{\rho }_{\mathrm{D}}^{\mathrm{MD}}\right)}{\ln \left(1+\mathsf{\Delta }{\rho }_{\mathrm{D}}^{\mathrm{MD}}\right)}\right)}^2\right]\right),$$

(21)

where $\left\{\mathsf{\Delta }{\sigma }_{11}^{\mathrm{MD}},\mathsf{\Delta }{\sigma }_{22}^{\mathrm{MD}},\mathsf{\Delta }{\phi }^{\mathrm{MD}},\mathsf{\Delta }{\rho }_{\mathrm{D}}^{\mathrm{MD}}\right\}$ are the ranges of the corresponding values from the MD, $N=180$ is the number of different MD systems of porous aluminum and different deformation paths at the same test temperature, and $M=60$ is the number of comparison points along each MD curve, arranged in steps of 0.005 of engineering strain. The total number $N·M$ of comparison points is about 10,000. Probability distributions in the parameter space are plotted, and an actual (optimal) set of model parameters is selected in the area of maximal probability as discussed in Section 3.2. The used parameters are collected in

Table 1. Parameters of the plasticity and fracture model: $\gamma$ is the surface tension coefficient, ${\epsilon }_{\mathrm{m}}$ is the dislocation multiplication energy from the dislocation multiplication term, $\alpha$ is the hardening coefficient from the Taylor hardening law, $k_{\mathrm{n}}$ is the coefficient of the energy barrier of dislocation nucleation from the nucleation term, and $k_{\mathrm{a}}$ is the coefficient of annihilation of dislocations.

Parameter	Value
$\gamma$	1 J/m2
${\epsilon }_{\mathrm{m}}$	0.8 eV/b
$\alpha$	0.5
$k_{\mathrm{n}}$	0.1
$k_{\mathrm{a}}$	200

.

3. Results

3.1. Dislocation Nucleation in Nanoporous Aluminum

Figure 7, Figure 8 and Figure 9 present our MD data on plasticity incipience in nanoporous aluminum with different shapes and sizes of pores as well as in a perfect single crystal of aluminum. The threshold of dislocation nucleation is plotted in the strain space $\left\{E_{11},E_{22}\right\}$ in Figure 7a, Figure 8a and Figure 9a and in the coordinates {pressure in the solid phase, shear stress in the solid phase} in Figure 7b, Figure 8b and Figure 9b. The MD data represented by symbols are supplemented by the results of the trained ANN represented by solid lines. The ANN is used to plot these lines as follows: a number of beams are emitted from the undeformed state with a small angular step in the strain space $\left\{E_{11},E_{22}\right\}$. The nucleation strain distance function $Q$ is traced by means of ANN calculation for the deformed states along each beam until the value of $Q$ becomes positive, which signals the dislocation nucleation. One can conclude a close fit of the ANN-based lines and MD points in Figure 7, Figure 8 and Figure 9, which validates the application of the ANN-based description of the plasticity incipience.

Figure 7 shows the influence of the pore size using an example of spherical pores in comparison with the pore-free samples. Figure 8 examines the influence of raised temperature on the plasticity incipience in a perfect single crystal and in nanoporous aluminum with 8 nm spherical pores. Figure 9 compares nanoporous samples with different shapes of pores at the same porosity of 3.3% (the equivalent pore size is 8 nm). Figure 9 also includes the data for a larger MD system with 8 nm spherical pores and porosity of about 0.41%, which are added to separately study the effect of pore size and porosity.

3.2. Verification of Plasticity and Fracture Model and Identification of Parameters

The probability distribution in the parameter space is shown in Figure 10; it is presented by cross-sections over pairs of parameters plotted using 40,000 random sets of the model parameters $\mathsf{A}=\left\{ \gamma ,{\epsilon }_{\mathrm{m}},\alpha ,k_{\mathrm{n}},k_{\mathrm{a}}\right\}$. It should be mentioned that these cross-sectional plots accumulate all points along the axes of three other model parameters in one point; therefore, the actual position of the maximum in the five-parameter space can be dismissed from what appears in these plots. At the same time, Figure 10 provides the basic information about the obtained probability distribution: not a strict maximum, but any set of parameters that are in the area of high probability can be chosen as the optimal one. One can see in Figure 10a that the area of high probability forms a stripe along $\gamma$, which means that the model results are insensitive to the surface tension for the considered pore sizes. Nevertheless, the maximum probability is located at 1 J/m², and we select this value as shown in Figure 10a and Table 1; this value is not contradictive to a solid metal. The hardening coefficient $\alpha$ in Figure 10b can be also selected in a relatively wide range, from 0.5 to 1, and we select 0.5 as the more common one. In contrast, the dislocation multiplication energy ${\epsilon }_{\mathrm{m}}$ and the nucleation barrier $k_{\mathrm{n}}$ are well-defined in narrow intervals in Figure 10, which shows the high sensitivity of the model to these parameters. The model provides good results for the annihilation coefficient $k_{\mathrm{a}}$ of 200 or higher, and we select the value of 200.

Figure 11 compares the MD and the model predictions with the selected parameters using the example of axial stress ${\sigma }_{11}$. Three random deformation paths differently directed in the strain space $\left\{E_{11},E_{22}\right\}$ are considered for each porous sample, and different sets of these paths are used for different shapes of pores in order to explore more widely in this figure the complex deformation behavior examined in our approach. The correspondence in the elastic parts of curves and in the moments of plasticity incipience is related to the quality of training of the ANN and seems to be quite good. The following plastic behavior is also adequately described by the trained model of plasticity and fracture. Both the sharp increase and the decrease in ${\sigma }_{11}$ after the beginning of plastic deformation are caused by the convergence of longitudinal stress and transverse stress due to the plastic relaxation. After the transient processes that follow the plasticity incipience, the pore growth stabilizes stresses close to a certain level, which depends on the strain hardening. Strain hardening is determined by the hardening coefficient and the dislocation density, which, in turn, depend on the dislocation kinetics and the related model parameters.

4. Discussion

Because we consider both nanoporous and pore-free aluminum, the examined deformation paths in Figure 1 go far beyond the threshold of plasticity incipience in nanoporous aluminum, which additionally provides a number of data on the inelastic behavior governed by the dislocation activity and pore (damage) growth at tension. Bubble plots in Figure 1 show that the plasticity incipience leads to a fast relaxation of high negative pressure at all-around tension in sector $E_{11}>0, E_{22}>0$ of the strain space, while it mostly leads to the relaxation of high shear stress under the condition of lateral tension $E_{22}>0$ combined with axial compression $E_{11}<0$. The latter is obvious because there is no substantial negative pressure in the combined compressed–tensile states in the sector $E_{11}<0, E_{22}>0$, while the shear stress is large under such conditions. Thus, the examined range of deformation paths includes quite deferent inelastic behavior of nanoporous metal, which is a challenge for a model description, and a substantial advantage of the developed model is that it adequately describes these paths.

Figure 2, Figure 3 and Figure 4 show that the plasticity starts from the dislocation nucleation on the pore surface with subsequent emission into the surrounding material. The initial dislocation semi-loop is formed by leading Shockley partial dislocation, which correlates with what is observed in the case of pore-free solid FCC metals [10,41,42]. The partial dislocations that leave behind stacking faults are shown in Figure 2, Figure 3 and Figure 4 as disordered plane regions. Later on, one can observe a transition to the split perfect dislocations, which can be seen as stripes of stacking faults in Figure 2, Figure 3 and Figure 4. It is remarkable that under the considered complex loading, the dislocation nucleation can start either on the rounded or on the flat surfaces of the cylindrical pores (see Figure 3). In the previously examined [30] cases of uniaxial and hydrostatic loading, the dislocations always nucleate on the rounded surfaces, which is explained in [30] by the stress concentration in the vicinity of rounded surfaces in contrast to the flat ones.

Figure 7 shows that the presence of even nanoscale pores drastically reduces the threshold of dislocation nucleation. Variation of the pore size and porosity in the range 0.42–27% influences little the plasticity incipience plotted in the $\left\{E_{11},E_{22}\right\}$ space in Figure 7a, whereas the influence is much more pronounced in the stress space in Figure 7b. It is explained by the following: larger pores require lower shear stress in the surrounding material to initiate dislocation nucleation and emission. At large porosities, this lower level of shear stress in the solid phase is reached approximately at the same level of strains $\left\{E_{11},E_{22}\right\}$ as at small porosities. In particular, it leads to non-monotonic size dependence of the nucleation threshold in the $\left\{E_{11},E_{22}\right\}$ space (see Figure 7a), whereas the size dependence in the stress space remains monotonic (see Figure 7b).

Experimental investigation [80] on the laser-induced shock loading of thin aluminum films gives the spall (tensile) strength of about 6.1 ± 0.5 GPa, which does not contradict our results presented in Figure 7b and corresponds to the shear stress level of about 1 GPa.

The temperature rise from 300 K to 700 K greatly reduces nucleation strains in the case of pore-free solid aluminum as shown in Figure 8a, whereas it has a little effect on the nucleation strains in the case of nanoporous aluminum. These dependencies replotted in the stress space in Figure 8b show a more significant difference due to the temperature-induced reduction of the elastic modules. The lower temperature sensitivity of the strength of nanoporous aluminum can be explained as follows. Thermal fluctuations increasing with temperature and the temperature-induced decrease in the shear modulus are the predominant factors for pore-free aluminum. The stress concentration and the presence of surfaces with partially broken interatomic bounds are the predominant factors for nanoporous aluminum, and these factors do not substantially depend on temperature.

Figure 9 plotted for the same volume of pores shows that the pore shape subtly effects the nucleation threshold with the weakest case of cylindrical pores and the strongest case of cubic pores. This order of materials in terms of strength corresponds to the conclusions of paper [30] about an easier nucleation of dislocation on the rounded surfaces in comparison with the flat ones. The ANN correctly predicts the nucleation strains for all considered pore shapes (Figure 9a). Although one can observe in Figure 9b some deviation in stresses for cylindrical pores and spherical pores at stress states close to hydrostatic loading, this deviation is a shift of no more than 0.1–0.2 GPa (about 5%) along the pressure axis, which is acceptable. MD data for a large MD system (100 × 100 × 100 lattice parameters or about 40 × 40 × 40 nm³) with 8 nm spherical pores are also plotted in Figure 9. One can conclude that the nucleation threshold in the stress space (Figure 9b) coincides with the cases of basic and large MD systems with the same pores. It means that the nucleation stress depends on the pore size and not on the porosity, which differs by a factor of eight between the two MD systems. The threshold in the strain space in Figure 9a is slightly displaced due to the influence of porosity.

We further develop the mechanical model proposed in [30] by taking into account the ANN-described threshold of dislocation nucleation and tensor equation of state. The developed machine learning-based model of plasticity and fracture is applied to describe the deformation behavior of various porous samples along multiple deformation paths in comparison with the MD. The incorporation of the ANN greatly increases the precision of the description of the elastic stage and of the moment of plasticity incipience. The obtained adequate correspondence for the following inelastic stage additionally substantiates the kinematics of deformation proposed in [30] for various shape of pores.

5. Conclusions

We performed a large number of MD simulations of tension of pore-free solid and nanoporous aluminum and their generalization was made in the form of a machine learning-based approach. This approach included an ANN for approximation of both the tensor equation of state in the elastic domain and the function describing the dislocation nucleation threshold. It also included a description of the following inelastic evolution of material by means of a mechanical model of plasticity and fracture, which is trained according to the MD results by means of the Bayesian algorithm of parameter identification. Our research revealed specific features of plasticity incipience in nanoporous aluminum with different pore shapes under complex deformation conditions and proposed an automated method of construction of the models of inelastic behavior of such metallic systems. The future development of this work supports application of the established model to describe the experimental results of the dynamic spall fracture of metals. Another possible direction of further development is consideration of other shapes of pores, such as octahedral nanovoids observed in quenched and aged aluminum [81].