Crystal Plasticity Parameter Optimization in Cyclically Deformed Electrodeposited Copper – A Machine Learning Approach

1. Introduction

Constitutive models enable to predict the behaviour of materials subjected to external loadings from small [1] to large [2] length scales. In some most straightforward problems, the constitutive models are simple and their parameters are relatively easy to establish using standard experiments. For example, establishment of the Young’s modulus of an isotropic, linearly elastic material is unambiguous using just the result of uniaxial tension test. However, in more complex situations, the constitutive models tend to use multiple non-linear equations and many parameters in order to capture the material’s behaviour. Such a situation is present e. g. in crystal plasticity, that is the topic of this paper.

The crystal plasticity (CP) theory in its most standard version accounts for dislocation glide on strictly defined crystallographic planes. Therefore, the information about orientations of crystallites in a polycrystalline microstructure of a given material has to be provided. Besides taking into account these microstructural details, the hardening equations of crystal plasticity are typically similar to conventional plasticity and typically account for isotropic and kinematic hardening. Various optimization algorithms such as gradient optimization [3], Newton-Raphson algorithm [4], Levenberg-Marquardt method [5], Bayesian optimization [6], particle swarm optimization [7] and evolutionary algorithms [8,9,10,11,12,13,14,15,16,17,18,19] are used in order to establish the correct set of material parameters (see the introduction in [19] for a thorough discussion). However, all of them share the basic disadvantage: one has to repeat the optimization in the case of having new material or new test result.

In order to overcome this limitation, Wang et al. [20] developed a neural network based methodology for parameter optimization. The idea behind the approach is performing the training of a neural network (NN) once and then building the database of predicted outputs using the trained NN. The created database is used to predict the set of optimal parameters based on the provided input. Both the neural network and the database can be stored and thus the main advantage of the proposed approach is that one does not need to perform a new optimization to get the parameters of a different material. Similar idea was presented in [21]. The main difference wrt. [20] is that machine learning was used only to establish reasonable parameter bounds and then fine-tuning was performed using a genetic algorithm.

In this article we propose to perform the optimization for the case of elasto-viscoplastic self-consistent (EVPSC) crystal plasticity model [14]. The investigated case is that of electrodeposited copper subjected to cyclic elasto-plastic deformation. We demonstrate that the proposed method is robust enough in order to provide correct predictions even in the highly complex problem as the one considered.

We build on the idea of using ML for the optimization of constitutive model parameters as described in [20]. However, there are important differences wrt. to op. cit. First of all, we study the case of a two-scale elasto-viscoplastic model, while the single-scale viscoelastic model was used in [20]. Second, we apply the LSTM approach directly to cyclic plasticity model data, while deep neural networks applied to data preprocessed using singular value decomposition technique were used in op. cit. Finally, single loading-unloading was considered in [20], while up to 50 cycles are simulated in the current paper.

The article is structured as follows. After this introductory section, we describe the methodology (Section 2). This includes the description of te EVPSC model (SubSection 2.1) that was already thoroughly described in [14,22] and machine learning (ML) methodology (SubSection 2.2). These are followed by the optimization results using two approaches for feature extraction from stress-strain curves (Section 3). Discussion (Section 4) is related to the applicability of the proposed method to other non-obvious tasks. The paper ends with conclusions (Section 5).

2. Methodology

This section contains two subsections. The first one (Section 2.1) covers the description of the constitutive model together with the corresponding micromechanical scheme, which was used to provide input data for training the machine learning framework. The second one (Section 2.2) describes the machine learning approach applied for the purpose of parameter optimization in this paper.

2.1. The EVPSC Model

The sequential elasto-visco plastic self-consistent (SEVPSC) code1 was recently developed and validated for the case of electrodeposited copper films [14] and is based on the sequential linearization scheme by [22]. The scheme coincides with the additive tangent approach [23] for the case of the Mori-Tanaka scheme. The results for the case of the self-consistent case were presented in [24] and [14]. Details of the self-consistent formulation can be found in the aforementioned papers. Nevertheless, some equations essential for the understanding of the present paper are provided in the following.

The strain rate, stress rate and stress have to fulfill the conditions stating that the average response over the polycrystal is consistent with the corresponding macroscopic values:

$$\dot{\overline{\epsilon }}=\langle {\dot{\epsilon }}_g\rangle ,\dot{\overline{\sigma }}=\langle {\dot{\sigma }}_g\rangle ,\overline{\sigma }=\langle {\sigma }_g\rangle .$$

(1)

In order to account for the interaction of the grain-level and overall quantities, the additive tangent interaction law proposed by [25] or [26] is adopted:

$${\dot{\epsilon }}_g-\dot{\overline{\epsilon }}=-{\overline{\mathbf{M}}}_{*}^{\mathrm{e}}·({\dot{\sigma }}_g-\dot{\overline{\sigma }})-{\overline{\mathbf{M}}}_{*}^{\mathrm{v}}·({\sigma }_g-\overline{\sigma }),$$

(2)

where $\dot{\overline{\epsilon }}$, $\dot{\overline{\sigma }}$ and $\overline{\sigma }$ denote the overall strain rate, stress rate and stress. In Eq.(2), ${\overline{\mathbf{M}}}_{*}^{\mathrm{v}}$ and ${\overline{\mathbf{M}}}_{*}^{\mathrm{e}}$ are the fourth order inverse Hill tensors for a purely viscous problem and a purely elastic problem, respectively.

The parameters to be established with the proposed ML methodology concern the isotropic and kinematic hardening at the level of individual grain. In the following the formulation of the single grain (subscript g) constitutive model is provided. The total strain rate is calculated by adding elastic ${\dot{\epsilon }}_g^{\mathrm{e}}$ and viscoplastic strain rates (the current version of the model is formulated in the small-strain theory) ${\dot{\epsilon }}_g^{\mathrm{v}}$:

$${\dot{\epsilon }}_g={\dot{\epsilon }}_g^{\mathrm{e}}+{\dot{\epsilon }}_g^{\mathrm{v}}$$

(3)

In order to account for the elastic part, the anisotropic Hooke’s law is used:

$${\dot{\sigma }}_g={\mathbf{A}}_g^{\mathrm{e}}·{\dot{\epsilon }}_g^{\mathrm{e}}\mathrm{or}{\dot{\epsilon }}_g^{\mathrm{e}}={\mathbf{M}}_g^{\mathrm{e}}·{\dot{\sigma }}_g,$$

(4)

where ${\mathbf{A}}_g^{\mathrm{e}}$ is the fourth order elastic stiffness tensor. The associated compliance tensor is ${\mathbf{M}}_g^{\mathrm{e}}={\left({\mathbf{A}}_g^{\mathrm{e}}\right)}^{-1}$.

We use rate-dependent CP formulation for a constitutive description at the single grain level. The plastic strain rate tensor at this level is calculated as a sum of shears on slip systems:

$${\dot{\epsilon }}_g^{\mathrm{v}}=\sum _{r=1}^N{\dot{\gamma }}^{r,g}{\mathbf{P}}^{r,g},\mathrm{where}{\mathbf{P}}^{r,g}=\frac{1}{2}({\mathbf{m}}^{r,g}\otimes {\mathbf{n}}^{r,g}+{\mathbf{n}}^{r,g}\otimes {\mathbf{m}}^{r,g})$$

(5)

where ${\dot{\gamma }}^{r,g}$ is the slip rate, $\{{\mathbf{m}}^{r,g},{\mathbf{n}}^{r,g}\}$ are the slip direction and the slip plane unit normal respectively. Superscripts r and g denote respectively a given slip systems and grain. ${\mathbf{P}}^{r,g}$ is the symmetric Schmid tensor. The resolved shear stress (RSS) is a projection of the Cauchy stress tensor on the given slip system:

$${\tau }^{r,g}={\mathbf{P}}^{r,g}:{\sigma }_g.$$

(6)

The slip rate depends on the RSS as well as on the slip system backstress ${\chi }^{r,g}$ and the critical resolved shear stress (CRSS) ${\tau }_{\mathrm{c}}^{r,g}$ [27]:

$${\dot{\gamma }}^{r,g}={\dot{\gamma }}_0{\left|\frac{{\tau }^{r,g}-{\chi }^{r,g}}{{\tau }_{\mathrm{c}}^{r,g}}\right|}^n\mathrm{sign}({\tau }^{r,g}-{\chi }^{r,g}).$$

(7)

Note that in Eq. 7 both isotropic and kinematic hardening is included. The other quantities are the reference slip rate ${\dot{\gamma }}_0$ and the inverse of the strain rate sensitivity n. The tangent viscous compliance tensor is equal to a derivative of viscoplastic strain rate over the stress and thus equals:

$${\mathbf{M}}_g^{\mathrm{v}}=n{\dot{\gamma }}_0\sum _{r=1}^N{\left|\frac{{\tau }^{r,g}-{\chi }^{r,g}}{{\tau }_{\mathrm{c}}^{r,g}}\right|}^{n-1}\frac{\mathrm{sign}({\tau }^{r,g}-{\chi }^{r,g})}{{\tau }_{\mathrm{c}}^{r,g}}{\mathbf{P}}^{r,g}\otimes {\mathbf{P}}^{r,g}.$$

(8)

Evolution of the CRSS and backstress are governed by isotropic and kinematic hardening laws, respectively. The linear-exponential isotropic hardening evolution looks as follows:

$${\tau }_{\mathrm{c}}^{r,g}={\tau }_{\mathrm{c}}^0+h_1{\gamma }^g+\left({\tau }_{\mathrm{sat}}-{\tau }_{\mathrm{c}}^0\right)\left(1-e^{-h_{\mathrm{b}}{\gamma }^g}\right),$$

(9)

where ${\tau }_{\mathrm{c}}^0,{\tau }_{\mathrm{sat}},h_1,h_b$ are material parameters of the isotropic hardening and ${\gamma }^g={\int }_0^t{\dot{\gamma }}^{g}dt$ is the accumulated shear in the grain. The Ohno-Wang kinematic hardening law [28] reads:

$${\dot{\chi }}^{r,g}=h_{\mathrm{k}}{\dot{\gamma }}^{r,g}-h_{\mathrm{kb}}\left|{\dot{\gamma }}^{r,g}\right|{\left(\frac{{\chi }^{r,g}{h}_{\mathrm{kb}}}{h_{\mathrm{k}}}\right)}^{h_{\mathrm{m}}},$$

(10)

where $h_{\mathrm{k}}$, $h_{\mathrm{kb}}$ and $h_{\mathrm{m}}$ denote the inematic hardening parameters. We have thus 4 isotropic hardening and 3 kinematic hardening parameters which yields in total 7 hardening parameters subjected to optimization, cf. [14].

2.2. Machine Learning

Long Short-Term Memory (LSTM) networks are a type of Recurrent Neural Networks (RNN) [29]. RNNs generate output by sequentially applying the same weights to each element in the input and network’s internal state. LSTM model architecture was selected because of its proven performance in other sequential prediction tasks, such as speech recognition and low parameter count when compared to fully connected networks. Our implementation is a unidirectional LSTM as defined in Pytorch 2.0 [30] with a hidden state vector length of 2048. All three curves from the dataset are generated simultaneously by the network. Output is a tensor whose shape depends on the approach. Two ways of evaluating the fitness of a given parameter set have been tested. The first was exactly the same as in [14], i.e., the stress values at the points where loading changes its direction were compared, cf. Figure 1. This approach will be from now on termed App1. In this case the shape of the output tensor is 2-by-100. In the second approach (App2), the stress-strain curves for each example were replaced with 3 simpler curves: stress amplitude A, mean stress $\overline{\sigma }$, and slope of each cycle $\alpha$ in order to simplify the prediction task, cf. Figure 2. In this second case the shape of the output tensor is 3-by-50. Input to the network is the same vector of seven parameters repeated 100 or 50 times (depending on the approach) to match the desired output sequence length.

Both input and output data were normalized to improve training stability by eliminating large differences in the scale of data. The model was trained for 300 epochs with a batch size of 16 and an initial learning rate of 1e-3. The learning rate was decayed for each batch with cosine annealing schedule. Warm restarts were applied after the first 50 epochs and after doubling the previous learning rate decay period after that until the end of training [31]. Weight decay of 1e-5 was used for regularization. Adam optimization algorithm was applied to train the network with mean squared error as a loss function. The architecture of the network is shown in Figure 3.

3. Results

The method was applied as follows. First, the database of SEVPSC code results was created by generating 36370 loading-unloading curves. Each result was produced by running the SEVPSC code for the same loading scheme but with a different set of parameters. The parameters were generated by selecting every possible combination of parameters from a list of values. The list for each parameter was generated by selecting 5 equally spaced values from a set range (cf. Table 1). Then, part of the database was used for training the LSTM network, while the remaining part was moved to the test set (the ratio of training set size to the whole set size was equal to 60%). The MSE obtained in both approaches are shown in Figure 4.

Next step was thus using the trained LSTM network for parameter optimization. The approach consists in using the trained LSTM to generate predictions in points lying in the parameter space that are more densely distributed than the set used for training. For this 35.8 million results were generated. Then, for the reference case (obtained by running the SEVPSC code with random parameters), the closest result in the LSTM-generated database is searched by calculating the mean squared error (MSE) between reference case curves and generated examples and selecting the generated example with the lowest MSE. We have tested the method for 30 different parameter sets. The 6 representative parameter sets obtained using this technique are shown in Table 2. The corresponding stress-strain curves obtained using the optimized parameters compared against the reference ones are shown in Figure 5. The entire group of parameter sets is presented in the Supplementary material. We have divided the results into 6 categories:

• Very good or reasonable agreement of SS curves obtained using parameters optimized in both apporaches – Figure 5 a and supplementary Figure S.1,
• Disagreement in the first cycle and reasonable agreement of SS curves obtained using parameters optimized in both apporaches – Figure 5 b and supplementary Figure S.2,
• Reasonable agreement of SS curves obtained using parameters optimized in App 1 (lack of convergence for App 2 parameters) – Figure 5 c and supplementary Figure S.3,
• Reasonable agreement of SS curves obtained using parameters optimized in App 2 (lack of convergence for App 1 parameters) – Figure 5 d and supplementary Figure S.4,
• Striking disagreement or lack of convergence – Figure 5 e and supplementary Figure S.5,
• Lack of convergence for the optimized parameters in both approaches – Figure 5 f and supplementary Figure S.6.

Out of 30 cases examined, 13 cases show very good or reasonable agreement in both approaches (category 1). 3 cases show reasonable agreement in late cycles but disagreement in the first cycle (category 2). It can be thus stated that more than half of the cases the method worked correctly in the case of both approaches. Then, in 4 cases App 1 provided correct prediction while App 2 did not converge (category 3). The reverse situation (category 4) was present in another 4 cases. Therefore, in 24 cases (80%) at least one approach provided the reasonably correct set of parameters. Concerning category 5, in 1 case, both approaches provided parameter sets leading to striking disagreement in SS curves (Figure 5e and Figure S.5.b). In 1 another case (cf. supplementary Figure S.5.a) there was striking disagreement for App 2 and lack of convergence for App 1. Finally, in 4 cases, neither approach provided a parameter set that can converge in actual SEVPSC simulation (category 6).

Let us now take a closer look on the obtained results. In terms of the approach, it seems that both of them are balanced and we cannot argue that either App1 or App2 is better. Cases where the predicted SS curve agrees with the actual one (Category 1, cf. Figure 5a and Figure S. 1.) prove the potential of the method. The fact that almost half of cases show only sligt deviations of the shape of the curve in the first cycle is actually a little bit surprising as the information about the shape of the first cycle was absent in both methods. It is thus easy to understand why 3 other cases (Category 2, cf. Figure 5b and Figure S. 2.) agree with the reference curve in the cyclic deformation regime while having some deviation in the first cycle. Morover, these 3 cases can be treated as correct prediction since the LSTM had no information about the shape of the initial loading curve. In the rest of cases (Categories 3-6) there was either agreement or lack of convergence, except for the curve shown in Figure 5e. Cases where the simulations with optimized parameters fails to converge are not useful but also harmless. Thus, only one case out of 30 led to a misleading prediction.

4. Discussion

The trained network was used for parameter identification. But it could be also used as a surrogate model enabling to obtain the same overall information in far shorter time. Such an approach was already presented in a number of papers. e.g., in [32] the surrogate ML model was trained using data from discrete dislocation dynamics. The network was able to tell what will be the stress-strain curve based on the initial dislocation configuration. On the other hand, Deshpande, Lengiewicz and Bordas [33] trained their deep learning (DL) network based on large deformation finite element method (FEM) simulations. Such an approach made it possible to predict the deformed shape based on the initial specimen geometry and applied boundary conditions.

Many surrogate models were also constructed in the framework of crystal plasticity (CP). A surrogate model for the ViscoPlastic Self-Consistent (VPSC) model was built in [34]. The researchers developed both forward and backward ML configurations. The first enabled to predict final crystallographic texture and stress-strain curve based on the initial texture and hardening parameters. The second enabled to predict initial texture and hardening parameters based on the final texture and stress-strain curve. A neural network enabling to predict stress-strain curve and texture but this time trained using the crystal plasticity finite element method (CPFEM) was presented in [35]. CPFEM trained ML architecture was used to predict cyclic stress values in low cycle fatigue experiments in [36]. An interesting contribution was presented in [37] where the ML was trained on elastic-viscoplastic fast fourier transform (EVPFFT) model and the prediction concerned spatially resolved crystallographic orientations. A novelty in [38] was application of CP-trained ML to non-monotonic strain paths. In [39] it was demonstrated that using an encoder-decoder deep learning framework it is possible to regenerate stress-strain curves for a material subjected to several loading-unloading cycles.

5. Conclusions

The paper demonstrated the possibility to use a machine learning framework based on LSTM to efficiently optimize material parameters in a highly complex situation of modelling the cyclic deformation of copper films using the elasto-viscoplastic self-consistent model with isotropic and kinematic hardening. One can conclude that the method is able to provide a reasonably accurate guess of the parameter set that in most cases provides the correct cyclic stress-strain.

The big advantage of the proposed machine learning-based method over optimization algorithms is that in the case of having a different material subjected to the same loading conditions one does not need to perform any other simulations to get the optimal parameters. One could argue that similar result could be obtained using only a large enough database. Using such a database instead of the machine learning approach has however two disadvantages. First, in order to have accuracy comparable to the one attained here, one should perform about 10 thousand times more simulations (because the number of simulations used for training was 3670 and the number of data points generated by trained ML framework was 35.8 million. Second, in order to use the database for future purposes, one would have to store the whole database of results (which would be 35.8 million times 977 kB = 35 TB), while in our case it is enough to store the trained model (which would be just 50.8MB) and when necessary regenerate the database with arbitrarily large or small number of example inputs to achieve a desired size of database.

Of course, the developed approach is not able to provide the material parameters based on different loading conditions than were used as an input. Therefore, in such a case one would have to repeat the optimization. The method was demonstrated for a highly challenging problem of cyclic deformation simulated with the micromechanical model accounting for both isotropic and kinematic hardening. It is expected that the approach is directly applicable to other cases. In particular, one can naturally apply the developed methodology to a case of studying low and high cycle fatigue with plasticity and damage models, cf. e.g., [40].