For the investigation in both applications, quasi-2D representative volume elements (RVEs) were generated using the DMG, which couples a particle simulation method with a radical Voronoi tessellation algorithm (Boeff, 2016). In the first step, the target grain size distribution is determined via a log-normal distribution. Hence, the average grain diameter as well as the standard deviation are required. With respect to prescribed distribution parameters, the number and size of spheres are predefined, which mimic the targeted grain size distribution. In the second step, spheres are randomly distributed into a finite volume which is larger than the intended final RVE. This finite volume is then compressed, allowing spheres to move freely under a repulsive potential and to avoid their overlapping. In the third step, updated sphere positions and diameters of each sphere from selected time steps are then fed to a radical Voronoi tessellation algorithm from the open-source software Voro++ (Rycroft, 2009) to construct RVEs. The resulting grain size distribution of these RVEs is then compared to the targeted grain size distribution, and the RVE with the minimum difference is selected accordingly. It must be noted that the shape of the RVE, generated using DMG, is rugged to leave the grain intact and to improve the mesh quality. In the forth step, to create the RVE for the microstructure simulations, the geometry of the 2D RVE is extruded for 1% of a side length of RVE and meshed with eight-nodes-linear-brick elements (C3D8) by using CUBIT (Sandia National Laboratories, 2016).

In the final step, periodic boundary conditions, following an approach introduced by Smit et al. (1998), are applied to the RVE. Further details on the implementation are described in Kulosa et al. (2017). The basic idea of this approach is that opposite nodes are coupled such that their displacements are the same. The global boundary conditions and strain are imposed to the reference vertex points V₁, V₂, V₄, and H₁, which are located at the corners of the RVE. An example of an RVE generated by using the introduced method is illustrated in Figure 1. Furthermore, comparison of diameter distribution between defined seed spheres and constructed RVEs with an average grain size μ of 6.0 μm and a standard deviation σ of 1.0 μm, and an average grain size μ of 13.0 μm and a standard deviation σ of 1.0 μm are illustrated in Figures 1C,D, respectively. From the comparison, grain size distributions of both RVEs are in good agreement with the targeted size distributions. In the next section, the crystal plasticity-based material model is described.

The material behavior of the FE simulation is described by a phenomenologically based crystal plasticity model. To resolve the heterogeneous deformation resulting from abrupt changes in mechanical behavior across grain boundaries of the considered polycrystal and to consider size effects between small and large grains, a nonlocal crystal plasticity model proposed by Ma and Hartmaier (2014) is implemented. As the applied nonlocal crystal plasticity model is already described in Ma and Hartmaier (2014), only an overview of the formulation is given. For further details on the non-local flow rule, the reader is kindly referred to Ma and Hartmaier (2014). In the following, quantities written in bold letters refer to vectors (small letters) and matrices of second rank tensors (capital letters). From the kinematics of deformation, the total deformation gradient F can be multiplicatively decomposed into the elastic deformation gradient F^e and the plastic deformation gradient F^p,

The elastic deformation is calculated using the Hooke's law. The plastic deformation is characterized by the plastic velocity gradient L^p, which is a function of the plastic deformation gradient F^p and its rate as,

For this study, a crystallographic slip of dislocations is defined as the only mechanism for plastic deformation. Thus, L^p is taken as the sum of the shear rates of all slip systems,

Here, γ∙α is the plastic shear rate. M_α = d_α⊗n_α is the Schmid tensor for slip system α, which is defined by the slip direction d_α and the slip plane normal n_α. The symbol ⊗ denotes the dyadic product of two vectors resulting in a second rank tensor. The total number of slip systems is N.

With respect to the nonlocal crystal plasticity model proposed by Ma and Hartmaier (2014), the flow rule and the hardening law can be expressed as:

and,

where, γ∙0 is the reference shear rate, and p₁ is the inverse value of the strain rate sensitivity. Furthermore, h₀ is the reference hardening parameter, χ_αβ is the cross hardening matrix, which is assigned as 1.0 for coplanar slip systems and 1.4 otherwise,τ^sat is the saturation slip resistance, and p₂ is a fitting parameter. The initial value of the slip resistance τ^α is defined as τ^0, and sgn() is a mathematical function that extracts the sign of a real number. The resolved shear stress τ_α for each slip system can be calculated from the stress S_α in the intermediate configuration or the state involving only the plastic deformation gradient F^p as,

The flow rule in Equation (4) consists of two additional back stresses τ^αGNDk and τ^αGNDi describing the hardening contributions from geometrically necessary dislocations (GNDs) (Ma and Hartmaier, 2014). The nonlocal constitutive model, in this context, is derived from the concept of super GNDs densities and incorporates the plastic strain gradient. Within a continuum mechanical approach, it is not possible to define crystallographic GND based on the Nye tensor in a unique way. To capture the internal stresses resulting from GND, the concept of super dislocations is followed, which allows us to define the dislocation Burgers vectors and line directions uniquely (Ma and Hartmaier, 2014). This hardening from plastic strain gradients is split up into an isotropic hardening part τ^αGNDi and a kinematic hardening part τ^αGNDk.

The second rank dislocation density tensor G in the reference configuration is computed from the curl of F^p as introduced by Nye (1953),

where Δ_jkl is the third rank permutation tensor and “l” represents the derivative with respect to the cartesian coordinate “l”. It must be noted that in Equation (7) the dislocation density tensor is written in index notation (G = G_ij). Since a reconstruction of meaningful crystallographic dislocation populations in a unique way is impossible, a unique definition of super GNDs is obtained by projecting the dislocation density tensor to the global Cartesian coordinates of the system. As a result, the stress fields of the crystallographic GNDs can be described with a good accuracy (Ma and Hartmaier, 2014), and the GND density tensor can be segmented into nine independent parts ρ̄α by evaluating,

where d_α and t_α are permutations of the Cartesian unit vectors as determined in Ma and Hartmaier (2014), and b is the magnitude of the crystallographic Burgers vector. The super GND densities for α = 1, 2, 3 represent screw-type superdislocations, while the remaining 6 components represent edge-type superdislocations, which are vital for determining the internal stress fields as a consequence of the super GNDs.

The isotropic hardening for the dislocation slip contributed by these super GNDs can be expressed using a Taylor-type equation,

Here, c₁ is the Taylor hardening coefficient or a geometrical factor [38], and μ is the shear modulus. χαβGND is the cross hardening matrix between crystallographic mobile dislocations and super GNDs.

The long-range internal stresses, caused by GNDs in dislocation pile-ups, contribute to the kinematic hardening effect. This part is calculated by evaluating the second order gradient of F^p, which results in a super GND gradient ρ_α in the form,

By evaluating these gradients within a small volume of dimension L³, the internal stresses S~GND in the intermediate configuration caused by dislocation pile-ups at grain boundaries can be calculated as explained in Ma and Hartmaier (2014). Thus, the kinematic hardening can be given by:

For the FCC crystal structure, the dislocation slip on the common crystallographic 〈110〉{111} slip systems is considered. On the other hand, we only take the dislocation slip on the crystallographic 〈111〉{110} slip systems into account for the case of BCC crystal structure.

In the previous section, the material model for the microscopic FE simulations was described. For the ML algorithms, homogenized values (or global values) that describe the RVE are used. In the following, the global homogenized parameters have the superscript RVE. The homogenization procedure is different for the two applications presented here (cf. sections 4 and 5). For the prediction of the damage evolution, the global values are homogenized according to the Hill-Mandel condition (Hill, 1963, 1972) in section 2.3.1 and with respect to the stiffness reduction in section 2.3.2. It is necessary to use such volume average technique, because the damage needs to be calculated locally. For the prediction of microstructural features from the flow curve, macroscopic stress and strain tensors are calculated with respect to the approach of Nemat-Nasser (Nemat-Nasser, 1999). In this case, we only need to formulate a macroscopic stress and strain tensor in order to calculate the flow curve. Therefore, we use a much simpler and numerically more effective efficient approach described in section 2.3.3. In this section, methods to homogenize global values or macroscopic properties from microstructure simulations are described.

Following the work of Hastie et al. (2008), SVR is an extension of linear regression and used for non-linear problems. In the following, the theory of SVR is briefly described. A more detailed description can be found in Appendix 7.3. The main idea is to gain a function fitting the given data points so that all points lie within a (small) distance of ϵ to the function (see Figure 3A). In Figure 3A, a simple two-dimensional problem is shown, in which all data points are supposed to be described by a linear function. The green area is called the margin, and its width is equal to two times ϵ. To obtain the best fit, the main task is to minimize the margin, and for doing so to solve a convex optimization problem (cf. Smola and Schölkopf, 2004; Hastie et al., 2008). Furthermore, so-called slack variables ξ are introduced, for measuring the relative distance by which the target distance of ϵ is violated (cf. Figure 3A). The points far away from the margin are the so-called support vectors. In addition, a regularization or cost parameter C is specified. It balances the contradictory goals of a good fit vs. a simple model by weighting the penalty for the slack variables. Furthermore, outliers have more influence in shaping the predicted output. To enable the algorithm to develop complex non-linear functions, so-called Kernels are introduced (Ng, 2016). Kernels are customizable to the needs of the target domain, which gives the algorithm the advantage to be adaptable to many problems. With the kernel function, it is possible to map the input data into an enlarged feature space. Since this mapping is in general non-linear, kernels enable SVRs to represent highly non-linear functions. In this work the Gaussian radial basis function kernel:

is used (Müller and Guido, 2017). The parameter γ controls the width of the Gaussian kernel. The decision function is then no longer linear, but rather a flexible weighted sum of Gaussian kernels.

RFRs are a combination of multiple Decision Trees (DTs) or, more precisely in our case, regression trees. It is a prototypical ensemble method, which builds a highly accurate predictive model by combining many simple models (often referred to as weak learners). Each DT predicts an output, and their results are averaged. DTs are hierarchy-based models where the goal is to find the right answer by “asking as few if/else questions” as possible (Müller and Guido, 2017). For regression, nodes contain the distinction whether a value is below or above a threshold value. The main idea is to split the feature space into regions using recursive binary partitioning (cf. Hastie et al., 2008), so that every new data point can be assigned to one region. A visual example of a RFR is given in Figure 3B. Here, the single DT has a so-called depth of two. The tree depth is equal to the longest number of consecutive nodes in a tree. Every DT starts with a root (the top node), which contains the first question, e.g., whether a chosen feature of the data point is smaller or larger than a specific value. The nodes of the last layer of the tree are called leaves. Each leaf corresponds to one target value, i.e., a single value of the output domain. Each data point is assigned to exactly one leaf by following the decisions down the tree. If a leaf contains only data points that correspond to the same target value, the leaf is called pure. Using DTs with pure leaves results in a model that can fit the training data perfectly, but can result in over-fitting. There are four important algorithm parameters that are tuned for the RFR in this work (cf. Müller and Guido, 2017). The number of used DTs (estimators) influences the amount of over-fitting and also the computation time. In addition, the maximum depth of each tree can be chosen specifically, or the tree is built until each leaf is pure or reaches a minimum number of samples inside the node. Furthermore, a criterion to decide whether to split a node needs to be defined, e.g., mean squared error. Another important parameter is the maximum number of features used for splitting a node. In general, a low value of this parameter means that each tree is different and may not need to be deep enough to be sufficiently accurate. A high maximum feature parameter or setting the value equal to the total number of features, results in DTs that are quite similar and thus defeating the purpose of an ensemble in the first place. The training data are fitted well by building deep trees and using the most distinctive features.

A variety of loading states are simulated to make the data base valid for damage occurring under general monotonous loading paths. Hence, nine displacement-controlled simulations with different loading states are performed: uniaxial tension, biaxial tension cases, and shearing as shown in Figure 1B. The nine loading cases are uniaxial tension in x- (Tx) and y-direction (Ty), biaxial tension Txy, T2xy, and Tx2y (see Figure 1B). In addition, four shearing cases were applied: Shearing in x- (Sx, S2x) and y-direction (Sy, S2y) according to Figure 1B. The RVE used for the creation of the data set for ML is presented in Figure 1A in section 2.1. It contains 51 grains with a mean grain size of 59μm and a standard deviation of 10μm, which results in a grain size range between 40 and 90μm. The material model, as well as the damage model and the homogenization methods, are described in sections 2.2 and 2.3, respectively. First, the local results are presented. Then the global material behavior is presented. As an example, Figure 4 shows contour plots of the von Mises stress and the damage parameter for uniaxial tension in x-direction.

In both Figures 4A,B, a strain localization in form of a band can be seen inside the RVE. Note that the damaged zone is split up because of the periodic boundary conditions. It is due to such a morphology of the damage band that a new homogenization scheme is required to homogenize it from the micro to the macro-scale (see section 2.3.2). At the flanks of the localization band, the stress is close to zero and the damage parameter has reached its maximum of 0.999. Furthermore, it is noted that the damage band propagates through the grains, i.e., in a transgranular manner, as one would expect for a ductile material, where damage and fracture occurs by void nucleation, coalescence and growth. From the simulations, relevant parameters for the homogenization are extracted: equivalent total, elastic and plastic strain, equivalent plastic strain rate, von Mises and hydrostatic stress, as well as the element volume. Locally, the parameters are computed as follows: The equivalent plastic strain is computed as described in section 2.2.1 and its rate is computed equivalently to the rate of the damage parameter according to Equation (20). The equivalent total and elastic strains are computed in the same way as the equivalent plastic strain using the Frobenius norm and the Green-Lagrange strain (cf. Equation 13). The total deformation gradient F is calculated as the gradient of the displacement, and the elastic deformation tensor is computed as F^e = FF^p−1 (Haupt, 2002). The von Mises and hydrostatic stress are computed according to Gross et al. (2011). The extracted values are homogenized as described in section 2.3.1 and 2.3.2: The global stress and strain values are the volume average of the local (element) values, and the global damage is calculated based on the stiffness reduction of the entire RVE. This results in eight global parameters: equivalent plastic strain (p) and its rate (p˙), equivalent total (ϵteq) and equivalent elastic strain (ϵeeq), von Mises stress (σ_vM), hydrostatic stress (σ_hyd), and the damage parameter (D) and its rate (D˙). After the homogenization, a further pre-processing of the points is applied (see Appendix 7.2), which spaces the data equally with respect to the equivalent plastic strain. Each data point represents one time step of the FE simulation. For each time step, there is a set of parameters consisting of the global parameters previously mentioned. Therefore, the complete data set has the size 9 × (●) ×8, where (●) is the number of time increments for each of the 9 loading cases applied to the single RVE (cf. Figure 1A), and 8 is the number of global parameters (p, p˙, σ_vM, σ_hyd, ϵteq, ϵeeq, D, D˙). In total, the time increments of all loading cases equal 3454. For all data points used in this application, the reader is kindly referred to Data Sheet 1_v1 in the Supplementary Material. The global values are used as the data set for the training and testing of the ML algorithms. The number of training and test data has been verified to be sufficient by using so-called learning curves, which are further described in section 4.2. The global material response in terms of the von Mises stress and damage rate with respect to the equivalent total strain can be seen in Figures 5A,B. The global behavior is given in the following by showing five out of the nine loading cases with the most significant difference in the material response.

It can be seen in Figures 5A,B that different loading conditions result in (quantitatively) different stress and damage evolution, although the general curve shapes are (qualitatively) similar. Each loading condition shows a distinct starting point for the initiation of damage, which corresponds to the maximum stress occurring at different global strains. In addition to the given global plots, it is worth having a look at the maximum global damage parameter. For the uniaxial and biaxial tension, the value is quite similar: 35.77% of maximum global damage for biaxial tension (Txy), and 35.78 and 37.7% for uniaxial tension in x-direction (Tx) and y-direction (Ty), respectively. The two shearing cases Sx and S2y have a lower maximum value for the global damage parameter. For shearing in x-direction (Sx), the maximum global damage occurring is 24.86%, and for S2y it is 27.18%. It should be noted that even though the tension cases share a similar maximum global damage value, the evolution of damage, with respect to the total equivalent strain, is different as seen in Figure 5B. In the following, the ML models and their parameters are presented.

For this application, Support Vector regression (SVR) and Random Forest regression (RFR) are used to predict the damage evolution D˙. The ML is conducted using python and scikit-learn 0.19.1 (cf. Pedregosa et al., 2011). The data set is split into training (75%) and testing (25%) data sets. For validation, the training set with 75% of the data is used as the “complete” data set, and therefore further split into a training set (for validation purposes) with 56.25% of all data points, and a validation set with 18.75% of the data. The testing and validation data set is unseen data that is only used for evaluating the final model, i.e., after the final ML parameters are set. The validation set acts as a test set during the fitting of the ML parameters. Before splitting the data into sets, the order of the data was randomized in a way that can be reproduced (constant random state of 666 Müller and Guido, 2017, Pedregosa et al., 2011). To assess the accuracy of the learning algorithms, the so-called R² score is used (Müller and Guido, 2017), which is computed as a fraction of the mean squared error and the variance (Pedregosa et al., 2011),

Here, y represents the output vector, index true indicates the reference output data, true, mean the mean value of the reference output data, and index pred represents the output data predicted by the ML algorithm. The total number of data points is assessed to be sufficient by using learning curves (cf. Ng, 2016) and cross-validation. Learning curves are a tool to check whether the number of data points used for training and testing is sufficient (Pedregosa et al., 2011). The training data is split several times into different set sizes to see the development in training and validation score with respect to the number of data used. Here, RFR is used to train and validate the model since its training process is very robust and shows only little sensitivity to the training parameters. As mentioned in section 4.1, a total number of 3454 data points are available. For the training and validation, 56.25 and 18.75% of the data is used, i.e., 2071 (training) and 519 (validation) data points, respectively. The training data is split seven times, so that the following absolute training split sizes result: 207, 517, 828, 1139, 1449, 1760, 2071. The validation set is 20% of each split. The resulting training and validation scores converge after using 1449 training data points to 97.3–97.5% for training and 82.3–82.8% for validation. Selecting the most predictive subset of features can help to avoid over-fitting. Therefore, the features were chosen according to conducted feature importance methods and to an ductile damage model from the literature (cf. Equation 25), which is formulated in a mathematically closed form as analytical function. Feature importance is used to assess the influence of each feature with respect to the result. The attribute importance can be understood as a value of how informative each feature is and therefore shapes the result. For the feature importance, RFR is used (cf. section 3.2) with the only non-default parameter being the number of Decision Tree (DT) (=500). Note that feature importance gives a rank of all features with respect to their impact on the results. Less important features are not necessarily trivial, and neglecting them does not automatically improve the results. Nevertheless, feature importance can provide an understanding of the relationship between input and output parameters with respect to the ML algorithms. As mentioned above, the training data contain 56.25% of the data, and the validation set accounts for 18.75% of all data points. The other 25% of the data is the test set, i.e., the unseen data that is only used after the training process to assess the ability of the ML algorithm to generalize. As mentioned in section 4.1, for the feature importance all extracted features are used without additional polynomial features or interactions. The results of the feature importance are presented in a bar plot in Figure 6.

The conducted feature importance results in the damage parameter being the most informative feature with a validation score of 89.37%. This leaves the importance of all other features to around 10% in total with 3.4% for the hydrostatic stress being the second most relevant feature, and the plastic strain rate the least important feature with 0.9%. Selecting half of the most important features (D, σ_vM, σ_hyd, ϵteq) produces a validation score of about 94.58%. One can see that the damage parameter and the stresses seem to be the most relevant, based on feature importance. For the RFR, selecting only the most relevant features (D, σ_vM, σ_hyd, ϵteq) leads to better results than other feature combinations. In contrast, these input parameters induce a lower accuracy for the SVR. Choosing the same features for SVR as chosen for RFR, results in a training score of just about 77.6% and a test score of 82.2%. Both score values are below acceptance. The SVR cannot extract enough information from the given features to approximate the damage rate sufficiently. Therefore, leaving out features causes an under-fitting problem so that all features are used: D, p, p˙, σ_vM, σ_hyd, ϵteq, ϵeeq. Taking a look at the analytical ductile damage model,

one can see a similarity in the input parameters compared to the results of the feature importance (Chaboche, 1988; Ambroziak, 2007). In the above Equation (25), s and S are material damage parameters. In addition, constant material parameters such as the Young's modulus [E = 228.96(GPa)] and Poisson's ratio [ν = 0.27(−)] are used, which can be calculated from uniaxial stress and strain curves. The input parameters for the analytical damage model are similar to the selected parameters by the feature importance, damage parameter and stresses. Later, the ML results are compared to the analytical damage model given in Equation (25) to investigate whether ML algorithms can describe the damage evolution at least as well as a well-established closed-form damage model.

Furthermore, cross-validation is used to find the best ML parameters. First, the most appropriate method to scale the data is determined for the SVR (as RFR does not require a scaling of the data). The input data are scaled according to a Gaussian normal-distribution with zero mean value and a variance of one (standard scaler Pedregosa et al., 2011). Moreover, cross-validation is used to assess the most suitable kernel and whether to use additional polynomial features. In this case, the Gaussian kernel and no additional polynomial features result in the highest accuracy. Furthermore, grid-search, i.e., finding a parameter set that results in the highest accuracy, is used to find the best parameter value of the regularization parameter C (cf. Equation 25) and the Gaussian kernel coefficient γ (cf. Equation 32). During grid-search, both parameters are fitted simultaneously. Here, the epsilon-SVR model is used, which is named after the parameter ϵ which can be found in Equation (28) of the Appendix 7.3. This precision parameter defines the distance between data point and target value, which is still considered accurate, and has no negative influence on the overall accuracy (Pedregosa et al., 2011). For the RFR, the cross-validation is used to choose the best number of DTs, the maximum tree depth and the split criterion of a node. A number of 500 DTs gives the best results with respect to a reasonable compromise on the computation time. Each DT is built until all leaves are pure, i.e., each last node corresponds only to a single target value, and the criterion to split a node is the mean absolute error. In Table 1, the optimum parameters of both ML algorithms are summarized. The other parameters, as defined in the scikit-learn library (Pedregosa et al., 2011), are set to their default values. RFR is rather robust with respect to the parameter values. Generally, SVR is more sensitive to parameter tuning. Therefore, its parameters were tuned within a smaller range. Within this range, the SVR parameters are not as sensitive to tuning. For example, changing the values for the parameters C and γ from their optimized values (based on the Grid-Search method) by 10% changes the training score by about 0.02% and the test score by around 0.04%. With the described data set and the fitted ML parameters, the two algorithms are trained. The results of the training processes are given in the next section 4.3.

For the homogenization of damage, two algorithms are used: SVR and RFR. The same randomly partitioned data set for the training and the testing process is used for both algorithms. The final training processes are conducted by using the previously defined parameters (cf. Table 1) and by using the features that lead to the best results as described in section 4.2. For SVR the training and testing processes both have a considerably high accuracy: 99.73% (training) and 98.25% (testing). One can see that the high accuracy for training as well as testing indicates no over- or under-fitting (high bias or high variance) problems. The same applies for the RFR: the training score is 97.66%, and the test score is 97.48%. The results of the algorithms are presented in Figures 7A,B. In both cases, only the testing data set (25% of all data) is displayed in the form of a predicted data against the target damage rate plot. The red lines in both figures represent the 5% mismatch area calculated based on the R² score (see Equation 24). Here, the predicted data are the output of the ML algorithm, and the real data are the reference damage rate gained from the simulations.

From the data set, one can see that the majority of data points have damage rate values below 1/s. Hence, the damage evolution is predicted more accurately for such values. Even though, both algorithms have a sufficiently high accuracy on the test set, the SVR is able to approximate the damage rate more accurately for higher damage rates. SVR lacks to approximate the damage rate sufficiently for values near zero as one can see in Figure 7A, where a reference damage rate of around 0.3/s is predicted for one point. The RFR has a lower accuracy for larger damage rate values but can approximate values near zero more accurately than the SVR. Some data points were predicted incorrectly with an error of more than 5% for both algorithms, but the SVR shows less scattering inside the ±5% mismatch area. Furthermore, the SVR is able to predict high values for the damage rate more precisely than the RFR.

In Figure 8A, the ML algorithm results are given in a damage rate against equivalent total strain plot for five loading cases. Both ML algorithms are able to capture the damage evolution with increasing strain, even though not every value can be predicted perfectly. Hence, the ML algorithms are capable of predicting the damage evolution for different loading states precisely. In general, the ML algorithms can approximate the material response with respect to damage behavior almost as well as the full-field FE simulations as shown in Figures 7, 8. The comparison of the trained ML algorithms to the analytical damage model is given in Figure 8B for the test data points for the uniaxial tension in x-direction. For the analytical damage model, the two material parameters had to be adjusted: s = 5.06(−) and S = 0.24(MPa) (cf. Equation 25). Both, SVR and RFR, are trained as described previously and shown in Figure 7. The ML algorithms are able to describe the damage evolution well as mentioned above. Nevertheless, SVR shows a slight over-fitting problem as the damage rate marginally decreases after the maximum of around 1.58(1/s) (cf. Figure 8B). One can see a small roughness in the course, but no over-fitting is visible for the RFR. Consequently, the RFR method is more robust to describe the damage evolution for the presented cases than the SVR method. Furthermore, the fitting of the algorithm parameter is less demanding for RFR compared to SVR. The analytical damage model is able to describe the general damage evolution (see blue line in Figure 8B). Nonetheless, some limitations for the analytical model are worth noting. According to the analytical model, a damage evolution is visible even before the actual damage initiation occurs (after about 10% of total strain). The reason for this is that in the numerical model we explicitly gave the limit of the strain value as the initiation criteria, while the damage evolves in the analytical model as soon as plasticity occurs; however, due to the selection of parameters it stays small until some level of plastic strain is reached. After the actual damage initiation, the analytical damage model also shows a gradual increase, although, some difference is observed between analytical model and the numerical simulations regarding the point of sharp increase in damage. Moreover, the analytical model is compared to only one loading state, and its generalization to a variety of loading states would require re-adjusting its material parameters. Furthermore, the analytical model does not allow to take microstructural quantities into account.

In a first step, 74 RVEs consisting of 100 grains with various grain size distribution parameters following a log-normal distribution function were generated using the dynamic microstructure generator (DMG) introduced in section 2.1. In this context, the average grain size μ and the standard deviation σ are varied between 6–13 and 0.1–1 μm respectively. To exclude any influence of crystallographic orientation on the deformation behavior of RVEs, 100 different sets of randomly chosen Euler angles have been assigned to all RVEs. In this way, the remaining factor influencing the strain hardening behavior must be the grain size distribution parameters of the microstructure. In the next step, the nonlocal crystal plasticity model described in section 2.2 is implemented onto a user-defined material model (UMAT) and applied in a finite element (FE) simulation with the commercial software ABAQUS to assess the mechanical response of RVEs. By using a nonlocal crystal plasticity model, size effects including the influence of grain size are taken into account. For this part of the study, a BCC crystal structure is assigned to all grains in the RVEs; nonlocal crystal plasticity parameters are summarized in Table 4 in the Appendix 7.1 (Vajragupta et al., 2017).

In the next step, the mechanical response of RVEs is simulated under a uniaxial tension loading condition, and macroscopic flow curves are homogenized from reference nodes using the method introduced in section 2.3.3. Examples of two RVEs with different grain size distribution parameters and the corresponding homogenized flow curves are shown in Figure 9. With a nonlocal crystal plasticity model, the influence of the grain size on the strain hardening behavior can be observed. These results prove the validity of implemented strain gradient crystal plasticity model and demonstrate that grain size effects can be incorporated properly in microstructure simulations. For the sake of simplicity, these flow curves are fitted with an empirical isotropic hardening law in order to reduce the dimensionality of the training data. In this context, the modified Voce law (Kim et al., 2013) is chosen and expressed as,

Y₀, R₀, R_inf, and β are material parameters to be determined, and p is the equivalent plastic strain. To parameterize the aforementioned hardening law from results of RVEs simulations, the nonlinear least square fitting method is implemented (Bates and Watts, 1988). As a result, two sets of calibrated modified Voce isotropic hardening parameters from two selected RVEs simulations are used to plot flow curves as illustrated in Figure 9. From the comparison, both fitted flow curves are in a good agreement with simulation results and can be used to represent microstructure simulations. Furthermore, the evolution of these fitted material parameters with respect to the average grain size is plotted as shown in Figure 10.

From Figure 10, the influence of the average grain size on fitted material parameters of the modified Voce law is observed. According to Equation (26), Y₀ is directly related to the yield stress. Fitted Y₀ as plotted in Figure 10A linearly decreases with an increasing average grain size, and standard deviation influences a scatter of Y₀ at the same average grain size. From Figures 10B,C, R₀ and R_inf non-linearly decrease with larger average grain size. These two parameters behave similarly to the Hall-Petch relation. However, the standard deviation does not contribute to a scatter of R₀ and R_inf. β, which inversely governs the slope of the hardening law and increases with an increasing average grain size. With respect to the hardening law, smaller average grain size results in a more pronounced strain hardening behavior. In the next step, these microstructure simulation results are fed as training data for ML models.

For this application, Support Vector regression (SVR) and Random Forest regression (RFR) are implemented to predict the average grain size producing a given material behavior, which is described by the parameters of the modified Voce hardening law. SVR and RFR are performed using Python and scikit-learn 0.19.1 (Pedregosa et al., 2011). The data are split into training (80 %) and testing (20 %) data sets. Similar to section 4, the R² score is used to evaluate the performance of ML models. To determine hyperparameters of selected ML models yielding the highest accuracy, Grid-Search with 3-fold cross validation is applied, which manually considers all combinations of hyperparameters in a search space.

In Table 2, the optimized parameters of both ML models are summarized while other parameters as introduced in the scikit-learn library (Pedregosa et al., 2011) are set to default values.

The training processes of both ML models for predicting the grain size from the flow curve are performed by using the defined parameters (cf. Table 2). For SVR, both, training and testing processes, give a high accuracy of 99.39% (training) and 97.95% (testing), respectively. These results indicate no over- or under-fitting issues. Similarly, trained RFR also results in a great accuracy for both training (99.62%) and testing (97.86%). The results of algorithms for the test data set (20% of all data) in the form of predicted grain size vs. reference grain size are shown in Figure 11. In this case, the predicted grain size data are the output from ML models and the reference grain size data are the grain sizes of RVEs used in microstructure simulations. From Figure 11, most of the data points from both trained ML models are within 5% error and there are only some data points, which give more than 5% error. Therefore, it can be concluded that there is no significant difference between both models in terms of scatter from the 100% accuracy line.

Furthermore, trained ML models are tested with data that are out of range of the training data. In this context, an RVE consisting of 100 grains with average grain size of 15 μm and a standard deviation of 0.1 μm are generated using DMG. Plastic behavior is again described with a nonlocal crystal plasticity model, with parameters given in Table 4 in the Appendix 7.1, and uniaxial tension loading conditions are simulated. This microstructure simulation is homogenized to obtain macroscopic flow curves and modified Voce law parameters are determined using the non-linear least square fitting method accordingly. The fitted modified Voce law parameters are summarized in Table 3. For the validation process, these parameters are then used as input for both trained ML models to determine the average grain size.

By comparing the average grain size of the RVE to produce the flow curve with grain sizes predicted from ML models, significant deviations are observed when out-of-range-data are used, because predicted grain sizes are always within the range of training data. Therefore, such results show that an application of these trained ML models are only valid when the input data are within a certain range. Furthermore, it must be verified that the output data lie within the space covered by the training data. To further improve accuracy and to extend applicability of trained ML models, more training data covering a wider range of grain sizes should be used. In any case, within the range of training data, predicted grain sizes are still in a very good agreement with the reference data.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.