A machine learning constitutive model for plasticity and strain hardening of polycrystalline metals based on data from micromechanical simulations

The deformation gradient is a fundamental concept in continuum mechanics, which describes how a material deforms or changes its shape under external forces or displacements. It quantifies the local deformation of a material at each point in space. The deformation gradient ${\boldsymbol{F}}$ is a tensor that relates the initial configuration of a material (reference state) to its deformed configuration (current state). As shown in figure 1, when considering an intermediate or relaxed state where all the material points are unloaded and only plastically deformed, the deformation gradient can be decomposed as:

$$\begin{align}{\boldsymbol{F}} = {\text{ }}{{\boldsymbol{F}}^{\,\text{e}}}{{\boldsymbol{F}}^p}\end{align} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

where ${{\boldsymbol{F}}^{\,\text{e}}}$ arises from the stretching and rigid rotation of the crystal lattice and ${{\boldsymbol{F}}^{\,p}}$ describes the plastic shearing on crystallographic slip systems [25, 26]. At low temperatures, it is assumed that the dislocation slip on well-defined crystallographic slip systems is the dominant source of plastic deformation; thus, according to Rice [27], the plastic flow can be described using the velocity gradient:

$$\begin{equation}{\boldsymbol{L}_p} = \left( {\boldsymbol{F}^p} \right)^{-1} \dot{\boldsymbol{F}^p} = \sum_{\alpha = 1}^N \dot{\gamma}^\alpha \left( {\boldsymbol{m}_0^\alpha} \otimes {\boldsymbol{n}_{0}^{\alpha}} \right)\end{equation} \tag{ 2 }$$

where ${\dot \gamma ^\alpha }$ is the shear rate of slip system $\alpha$, $N$ is the number of active slip systems, ${\boldsymbol{m}}_0^\alpha$ and ${\boldsymbol{n}}_0^\alpha$ are the slip direction and normal to the slip plane, respectively. As the crystal deforms, vectors connecting lattice sites are stretched and rotated according to ${{\boldsymbol{F}}^{\,\text{e}}}$, the slipping direction and normal to slip plane in the current deformed configuration can be characterized as:

$$\begin{align} \begin{array}{*{20}{l}} {\boldsymbol{m}^\alpha = \boldsymbol{F}^{\,\text{e}} \boldsymbol{m}_0^\alpha} \\ {\boldsymbol{n}^\alpha = \boldsymbol{n}_{0}^{\alpha} (\boldsymbol{F}^{\,\text{e}})^{-1}} \end{array} \end{align} \tag{ 3 }$$

The resolved shear stress on slip system $\alpha$ with respect to the reference configuration can be written as:

$$\begin{align}{\tau ^\alpha } = {\boldsymbol{S}} \cdot \left( {{\boldsymbol{m}}_0^\alpha \otimes {\boldsymbol{n}}_0^\alpha } \right)\end{align} \tag{ 4 }$$

${\boldsymbol{S}}$ is the second Piola–Kirchhoff stress and is defined in the intermediate configuration. For face-centered cubic metallic crystals, the kinetic law of a slip system can be described as:

$$\begin{align}{\dot \gamma ^\alpha } = {\text{ }}{\dot \gamma _0}{\left| {\frac{{{\tau ^\alpha }}}{{\tau _c^\alpha }}} \right|^{\frac{1}{m}}}{\text{sgn}}\left( {{\tau ^\alpha }} \right)\end{align} \tag{ 5 }$$

where ${\dot \gamma ^\alpha }$ is the shear rate on the slip system $\alpha ,$ which is subjected to resolved shear stress ${\tau ^\alpha }$ with a resistance of $\tau _c^\alpha$. ${\dot \gamma _0}$ is the reference shear rate, and $m$ is the slip rate sensitivity. Both are materials parameters [25, 26, 28].

The hardening behavior of slip system $\alpha$ is influenced by any slip system $\beta$ as follows:

$$\begin{align}\dot \tau _c^\alpha = {h_0}{q_{\alpha \beta }}{\left( {1 - \frac{{\tau _c^\beta }}{{{\tau _s}}}} \right)^p}\left| {{{\dot \gamma }^\beta }} \right|\end{align} \tag{ 6 }$$

where ${h_0}$ and $p$ are material parameters, ${\tau _s}$ is the saturated slip resistance due to dislocation density accumulation and ${q_{\alpha \beta }}$ is the cross-hardening matrix.

In this work, we employed the CPFEM implemented in the commercial software package ABAQUS [29] to run simulations capturing the stress–strain responses of polycrystalline materials under various monotonic loading conditions. The outcomes of these simulations provided the data necessary to construct the training dataset, which subsequently informed our ML model. The CPFEM implementation requires steps such as the creation of RVEs, the assignment of texture, and the careful definition of boundary conditions and application of the loads. These steps are outlined in detail in the following sections.

2.2.1. RVE generation

Classical continuum mechanics deals with the mechanical behavior of materials modeled as a continuum rather than discrete particles, hence neglecting the atomic structure of the materials. Additionally, it does not consider the fact that the mechanical properties can be size-dependent, altering significantly when considering extremely small scales. On the other hand, micromechanics explores the behavior of materials on the microscale, considering the influences of individual grains, phases, texture, or other microstructural components on the mechanical response. A fundamental concept in micromechanics is the RVE or the repeating unit cell. This element provides a micro-to-macro bridge, enabling us to represent the intricate details of the microstructure and to to operate within the simpler framework of continuum mechanics [30, 31].

In this study, we utilized the RVE concept to capture the complex behavior of a polycrystalline material. The RVE generated in ABAQUS consists of 343 distinct cells, with each cell representing an individual grain within the polycrystalline structure. This number of cells was specifically chosen to ensure that the RVE is large enough to accurately represent the complex behavior of polycrystalline materials, providing a realistic and comprehensive portrayal of their microstructural heterogeneity. Each individual cell has dimensions of 0.095 mm. Each of these cells was subdivided into 8 smaller elements, resulting in a total of 2744 elements in the RVE. The choice to subdivide each grain in the RVE into 8 smaller elements was made to capture detailed microstructural characteristics while maintaining computational efficiency. Each cell was assigned unique material properties, including crystallographic texture, to create a detailed microstructural representation of the material. The generated RVE is shown in figure 2.

Zoom In Zoom Out Reset image size

2.2.2. Texture

2.2.3. Periodic boundary conditions

Periodic boundary conditions (PBCs) are often employed in the simulations of physical phenomena within systems that are conceptually infinite or substantially larger than the scale of internal interactions. These boundary conditions play a crucial role in determining the realistic deformation behavior of RVEs. In PBCs, the displacements of opposing nodes on the RVE are coupled in such a way that the deformation can still represent an average strain over the RVE if demanded by the imposed deformation [35, 36]. To implement PBC, the approach introduced by Smit et al [37] was adopted, which involves constraining the representative pair of boundary nodes located on opposite faces. To achieve this, the equations proposed by Boeff [36] were utilized. As illustrated in figure 3, the RVE is characterized by eight nodes labeled as F1, F2, F3, F4, B1, B2, B3, and B4. Only F1, F2, B1, and F4 are independent nodes to which global boundary conditions are applied, as detailed in table 1. The remaining nodes F3, B2, B3, and B4, are dependent on the displacements of the independent nodes through coupling equations:

$$\begin{align}\begin{array}{*{20}{l}} {{u^{B3}} - {u^{F3}} = {u^{B1}} - {u^{F1}}} \\ {{u^{B2}} - {u^{F2}} = {u^{B1}} - {u^{F1}}} \\ {{u^{B4}} - {u^{F4}} = {u^{B1}} - {u^{F1}}} \\ {{u^{F3}} - {u^{F4}} = {u^{F2}} - {u^{F1}}} \end{array}\end{align} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

Table 1. The independent nodes within the RVE to which global boundary conditions are applied.

Nodes	X-direction	Y-direction	Z-direction
F1	Fixed	Fixed	Fixed
F2	Concentrated load	Fixed	Concentrated load
B1	Fixed	Concentrated load	Concentrated load
F4	Concentrated load	Concentrated load	Fixed

2.2.4. Loading condition

Elastic deformation is the reversible change in the shape and volume of a material under stress, which disappears once the stress is removed. In contrast, plastic deformation is a permanent shape change that remains even after stress is alleviated. The yield function serves as a quantitative threshold marking the transition from elastic to plastic deformation, which occurs when the applied stress reaches the yield strength of the material. The mathematical representation of this framework can be described as:

$$\begin{align}f\left( {\boldsymbol{\sigma}} \right) = {\sigma _{{\text{eq}}}}\left( {\boldsymbol{\sigma}} \right) - {\sigma _y} \end{align} \tag{ 8 }$$

In this formulation ${\sigma _{{\text{eq}}}}\left( {\boldsymbol{\sigma}} \right)$ is the equivalent stress function of the current stress tensor ${\boldsymbol{\sigma}}$ and ${\sigma _y}$ represents the yield stress of the material. When the yield function responds with a value equal to or larger than zero, it signifies that the material is undergoing plastic deformation. On the other hand, if the yield function returns a negative value, it indicates elastic deformation. Upon yielding, the plastic deformation in many metals is governed by the associated flow rule. This rule postulates that plastic strain develops in the direction of the gradient of the yield function, which is described as:

$$\begin{align}{\dot\varepsilon _p} = \dot \lambda \frac{{\partial f\left( {\boldsymbol{\sigma}} \right)}}{{\partial {\boldsymbol{\sigma}} }}\end{align} \tag{ 9 }$$

In this equation ${\dot \varepsilon_p}$ represents the plastic strain rate and $\dot \lambda$ is the plastic multiplier. The associated flow rule implies that the material will deform plastically in the direction where the yield function increases most rapidly. To maintain consistency with the yield condition, the plastic multiplier is derived such that the consistency condition $\frac{{{\text{d}}f}}{{{\text{d}}t}} = 0$ is satisfied during plastic loading. This ensures that, as the material deforms plastically, its stress state always remains on the yield surface, even if the yield surface itself evolves (e.g. due to hardening) [4, 38]. Oftentimes, a return mapping algorithm is applied to calculate the proper magnitude of the plastic strain increment in an iterative way.

As suggested by Hartmaier [23], the yield function can be perceived as a decision boundary within a classification problem. In machine learning, classifiers are designed to discriminate between two or more classes of data. Like a classifier, the purpose of the yield function is to separate two distinct regions in the stress–strain space. Drawing parallels between these two concepts, SVC can be employed as a powerful tool to represent the yield function. The primary objective is to precisely identify an optimal hypersurface that serves as the boundary between the elastic and plastic regions, analogous to the role of a yield function. SVC is known for handling high-dimensional datasets well. Its main strength lies in the convex nature of the optimization problem it solves, ensuring a consistent and unique solution and avoiding issues like local minima [39–41]. It is important to note that this optimization convexity does not automatically ensure the convexity of the decision boundary, which is paramount when modeling yield functions. A convex yield surface guarantees a unique and stable material response under any given loading condition. This convexity ensures physically consistent behavior, preventing any unpredictable transitions or non-physical manifestations [42]. The requirement for convexity in yield surfaces has been emphasized and derived in various works. Specifically, the work assumption presented by Naghdi and Trapp [43] and the assumption of maximum rate of dissipation by Rajagopal and Srinivasa [44] both underscore the importance and necessity of convexity in yield surfaces. SVC has a remarkable ability to maintain and reflect the convexity inherent in the training data when determining the decision boundary, particularly when trained on data with inherent physical convexity. Hence, SVC naturally captures and mirrors this characteristic without the need for additional convexity enforcements which makes it an ideal ML algorithm to train a model that can effectively serve as the yield function.

Shoghi and Hartmaier [24] demonstrated the functionality of an SVC-based yield function using 6-dimensional stress tensors as training data in the case of ideal plasticity. In the present work, as a significant advancement, we extended the feature space of the ML model to include both stress and strain components. The training data we employ now comprise the flow stresses and their corresponding plastic strain components. This results in a 12-dimensional feature vector, with six of these dimensions dedicated to stress and the remaining six to strain. In mechanics, six independent variables—comprising three normal and three shear components—are sufficient for describing the full stress or strain tensor. This enhancement broadens the capability of the model to capture strain hardening and deepen its understanding of material behavior. This ML yield function denoted as ${f_{{\text{ML}}}}\left( {{\boldsymbol{\sigma}} ,{\boldsymbol{\varepsilon}} } \right)$ is designed and trained such that, when presented with a stress and corresponding strain tensors, it can provide a response of either ${f_{{\text{ML}}}}\left( {{\boldsymbol{\sigma}} ,{\boldsymbol{\varepsilon}} } \right) = - 1$ or ${f_{{\text{ML}}}}\left( {{\boldsymbol{\sigma}} ,{\boldsymbol{\varepsilon}} } \right) = + 1$. These numerical responses indicate the elastic or plastic behavior of the material respectively. Beyond merely acting as a classifier, the ML yield function can also be employed as a basis for a data-oriented flow rule. Once yielding happens for every specific plastic strain, the function aims to identify a stress state where the yield function is zero, ensuring that the calculated flow stress aligns with the anticipated plastic deformation. This methodology is grounded in the consistency condition inherent in classical plasticity theory. The decision function, handling a 12-dimensional feature vector, is expressed as follows:

$$\begin{align}{f_{{\text{ML}}}}\left( {{\boldsymbol{\sigma}} ,{\boldsymbol{\varepsilon}} {\text{ }}} \right) = \mathop \sum \limits_{k = 1}^{{N_{{\text{SV}}}}} {\alpha _k}{y_k}\psi \left( {\left( {{\boldsymbol{\sigma}} ,{\boldsymbol{\varepsilon}} } \right),\left( {{{\boldsymbol{\sigma}} _k},{{\boldsymbol{\varepsilon}} _k}} \right)} \right) + b\end{align} \tag{ 10 }$$

${N_{{\text{SV}}}}{\text{ }}$ is the number of support vectors, which are a subset of the training data points that are crucial for defining the decision boundary. ${\alpha _k}{\text{ }}$ and ${y_k}{\text{ }}$ are coefficients associated with each support vector. ${\alpha _k}$ is the Lagrange multiplier or weight associated with the $k{\text{th}}$ support vector and ${y_k}{\text{ }}$ is the corresponding class label (+1 or −1). $\psi \left( {\left( {{\boldsymbol{\sigma}} ,{\boldsymbol{\varepsilon}} } \right),\left( {{{\boldsymbol{\sigma}} _{\boldsymbol{k}}},{{\boldsymbol{\varepsilon}} _{\boldsymbol{k}}}} \right)} \right)$ is the kernel function that measures the similarity or proximity between the input vector $\left( {{\boldsymbol{\sigma}} ,{\boldsymbol{\varepsilon}} {\text{ }}} \right)$ and each support vector $\left( {{{\boldsymbol{\sigma}} _{\boldsymbol{k}}},{{\boldsymbol{\varepsilon}} _{\boldsymbol{k}}}} \right)$. The kernel function is responsible for transforming the input space into a higher-dimensional feature space where the data are linearly separable. For the nonlinear problem at hand, the radial basis function (RBF) kernel $\psi \left( {\left( {{\boldsymbol{\sigma }},{\boldsymbol{\varepsilon }}} \right),\left( {{{\boldsymbol{\sigma }}_{\boldsymbol{k}}},{{\boldsymbol{\varepsilon }}_{\boldsymbol{k}}}} \right)} \right) = {\text{ exp}}( - \gamma \| ( {{\boldsymbol{\sigma }},{\boldsymbol{\varepsilon }}} )-$ ${\text{ }}( {{{\boldsymbol{\sigma }}_{\boldsymbol{k}}},{{\boldsymbol{\varepsilon }}_{\boldsymbol{k}}}} ) \|^2 )$ is chosen, and $b$ is the bias term, which is a constant offset that helps shift the decision boundary. The RBF kernel was selected due to its capability to effectively handle complex and non-linear patterns, especially in the context of our 12-dimensional fracture vector. In contrast to the polynomial kernel, which necessitates a predetermined degree of relationship, the RBF kernel offers flexibility and is often more proficient in such high-dimensional scenarios. The SVC model is trained through a dual optimization problem, which focus on finding the optimal values for the Lagrange multipliers ${\alpha _k}$. The dual optimization problem can be formulated as a maximum search of the expression

$$\begin{align}\mathop {{\text{max}}}\limits_\alpha \left( {\mathop \sum \limits_{k = 1}^{{N_{{\text{SV}}}}} {\alpha _k} - \frac{1}{2}\mathop \sum \limits_{k = 1}^{{N_{{\text{SV}}}}} \mathop \sum \limits_{l = 1}^{{N_{{\text{SV}}}}} {\alpha _k}{\alpha _l}{y_k}{y_l}\psi \left( {\left( {{{\boldsymbol{\sigma}} _k},{{\boldsymbol{\varepsilon}} _k}} \right),\left( {{{\boldsymbol{\sigma}} _l},{{\boldsymbol{\varepsilon}} _l}} \right)} \right)} \right)\end{align} \tag{ 11 }$$

with respect to the Lagragian multipliers under the constraints:

$$\begin{align}\begin{array}{*{20}{l}} {\mathop \sum \limits_{k = 1}^{{N_{{\text{SV}}}}} {\alpha _k}{y_k} = 0{\text{ }}} \\ {0 \le {\alpha _k} \le C} \end{array} \end{align} \tag{ 12 }$$

Here, $C$ is the regularization parameter limitting the magnitude of the Lagrange multipliers and, thus, controlling the trade-off between achieving a low training error and ensuring that the model generalizes well to unseen data. Once the dual optimization problem is solved, i.e. optimal ${\alpha _k}$ values are found, the decision function (10) can be used to make predictions [45, 46].

As the ML yield function is defined as a convolution sum over the support vectors, the gradient to the SVC decision function can be calculated as

$$\begin{align}\frac{{\partial {f_{{\,\,\text{ML}}}}\left( {{\boldsymbol{\sigma }},{\boldsymbol{\varepsilon }}{\text{ }}} \right)}}{{\partial {\boldsymbol{\sigma }}}} = \sum\limits_{k = 1}^{{N_{{\text{SV}}}}} - 2\gamma {\alpha _k}{y_k}\exp \left( { - \gamma {{\left\| {\left( {{\boldsymbol{\sigma }},{\boldsymbol{\varepsilon }}} \right) - \left( {{{\boldsymbol{\sigma }}_{\boldsymbol{k}}},{{\boldsymbol{\varepsilon }}_{\boldsymbol{k}}}} \right)} \right\|}^2}} \right)\left( {\left( {{\boldsymbol{\sigma }},{\boldsymbol{\varepsilon }}} \right) - \left( {{{\boldsymbol{\sigma }}_{\boldsymbol{k}}},{{\boldsymbol{\varepsilon }}_{\boldsymbol{k}}}} \right)} \right)\end{align} \tag{ 13 }$$

from which the plastic strain increments can be derived directly using a standard finite element formulation.

$C$ and $\gamma$ are two important hyperparameters in SVC with the RBF kernel. Parameter $C$ is the regularization parameter that controls the trade-off between maximizing the margin and minimizing the misclassification of the training examples. A lower $C{\text{ }}$ value allows for more misclassifications, leading to a broader margin, which can enhance generalization. However, a higher $C$ value imposes a stiffer penalty on misclassifications, resulting in a more rigid decision boundary characterized by a narrower margin. The parameter $\gamma$ defines the spread of the RBF kernel. It controls the extent of influence for each training instance, affecting the smoothness and flexibility of the decision boundary. A smaller $\gamma$ value results in a more localized influence, where each data point has a narrower impact range. This often leads to simpler and smoother decision boundaries. Conversely, a larger $\gamma$ value increases the range of influence for each training instance, resulting in a more complex decision boundary that conforms closely to the training data [40, 47]. The effects of different combinations of $C$ and $\gamma$ are listed in table 2. Each case in the table represents a combination of high or low $\gamma$ and $C$ values, showcasing the resulting characteristics of the decision boundary.

Table 2. The effects of different combinations of $\gamma$ and $C$ values on the decision boundary in SVC with the RBF kernel.

${\boldsymbol{\gamma}}$	${\boldsymbol{C}}$	Effect on Decision Boundary
Low	Low	A wider margin and a smoother decision boundary. It allows more misclassifications and tends to generalize better but might sacrifice accuracy.
Low	High	The margin becomes narrower, and the decision boundary becomes stricter. It penalizes misclassifications more severely, leading to a decision boundary that fits the training data more tightly.
High	Low	With a higher gamma value, the influence of each training example increases, resulting in a more complex and wiggly decision boundary. However, a low $C$ value still allows some misclassifications.
High	High	Results in a decision boundary that is both complex and strict. The decision boundary tightly fits the training data and minimizes misclassifications but might be wiggly.

It is important to note that the optimal choice of $C$ and $\gamma$ depends on the dataset and specific problem. Achieving the best performance often requires optimization techniques, such as grid search or Bayesian optimization, to determine the optimal hyperparameter values.

The resulting data from the CPFEM simulation serve as the basis for the generation of training data necessary to inform the ML model. However, before utilizing these data for training, a preprocessing step is necessary to classify the data into distinct classes representing the elastic and plastic regions of the stress–strain space. This feature extraction process can be carried out by utilizing flow stresses and their corresponding plastic strain components. For each equivalent plastic strain level, the flow stresses are upscaled towards the plastic region and downscaled to the elastic area. Following this methodology, the upscaled and downscaled flow stresses are then labeled, indicating whether they fall within the elastic or plastic class. These systematically labeled stresses, along with their corresponding plastic strain tensor, serve as the ground truth for our ML model.

During the preprocessing step, the critical task of upscaling and downscaling $N$ data points is performed. The size of the training dataset is directly determined by the value of $N$. As shown in figure 4 this operation occurs within the two ranges of [${F_e}$, ${C_e}$] and [${C_p}$, ${F_p}$]. In these ranges, ${C_e}$ and ${C_p}$ are the boundaries proximate to the actual flow stress in the elastic and plastic domains respectively. Meanwhile, ${F_e}{\text{ }}$ and ${F_p}$ define the outermost limits of these ranges. The real flow stress, marked in yellow, is set to a value of 1 in this schematic representation. The primary function of this real flow stress is to establish a convenient reference point, facilitating the subsequent generation of labeled training data in the elastic and plastic regions of the stress space, respectively. The correct configuration of the preprocessing parameters is essential for the accuracy and predictive capability of the ML model. Figure 5 illustrates the effects of varying preprocessing parameters on data transformation through three plots. Figure 5(a), set with a constant transformation range (${C_e}$, ${F_e}$, ${C_p}$, and ${F_p}$), illustrates the influence of varying $N$, the number of upscaled and downscaled points. This variability in $N$ directly influences the training size, which in turn can affect the overall quality of model training. Figure 5(b) shows the effect of varying ${F_e}$, and ${F_p}$ while $N$, ${C_e}$, and ${C_p}$ held constant. This visualization offers insights into the varied possibilities when the far boundaries in the elastic and plastic regions are adjusted, consequently altering the transformation range. Figure 5(c) depicts the potential transformations resulting from variations in ${C_e}$ and ${C_p}$ with $N$, ${F_e}$, and ${F_p}$ held constant. This illustration highlights the potential variations when adjustments are made to the outer boundaries of both the elastic and plastic regions, which in turn shifts the transformation range. It must be noted that for illustrative clarity, this visualization is simplified. In practice, the process is executed on 6-dimensional arrays of flow stresses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This visualization underscores the wide array of parameter combinations that can be explored during the preprocessing stage. By leveraging symmetry, a constant data generation range can be mirrored across both the elastic and plastic regions, thereby reducing the number of parameters that need determination. Utilizing techniques such as Design of Experiments (DoE) [48] facilitates systematic parameter exploration, aiming for optimal model performance. It is important to note that the preprocessing parameters also affect the hyperparameters $C$ and $\gamma$ of the SVC model. To identify this optimal set, it is essential to thoroughly explore the space of $C$ and $\gamma$ using grid search techniques. The grid search employs a cross-validation technique that helps to detect and mitigate overfitting by validating the model's performance on unseen data. Each $C$ and $\gamma$ combination is evaluated using a $k$-fold cross-validation process, where the dataset is divided into $k$ subsets. For each iteration, one subset serves as the test set, while the remaining $k - 1$ subsets are combined to form the training set [40, 49].

The feature vector powering the ML model is 12-dimensional. This vector comprises six stress components and six strain components, where the strain components represent only the plastic portion of the total strain. The emphasis on plastic strain is important because it represents how the material undergoes permanent deformation when subjected to stress. Such a focus is essential for accurately capturing and describing the strain-hardening behavior of the material.

In micromechanical simulations, the plastic strain starts accumulating right from the onset of deformation. However, it is a widely accepted convention to consider the transition from elastic to plastic behavior when the equivalent plastic strain reaches a small but non-zero threshold of 0.2%. This convention helps distinguish between the initial elastic deformation and the onset of more significant, irreversible plastic deformation. On the other hand, continuum mechanics adopts a different perspective. Within this framework, a distinct transition from elastic to plastic behavior is envisioned at a critical stress, below which the plastic strain is identically zero. This idealized yield onset, typically considered in continuum plasticity, enforces a clear transition between purely elastic behavior and the beginning of plastic deformation. In designing the ML model, alignment with the principles of continuum mechanics is intended, with zero plastic strain being assumed in the elastic regime up to the onset of plastic yielding. To reconcile this difference, the plastic strain data from the micromechanical simulations is adjusted to align with this assumption. As illustrated in figure 6, plastic strain data is adjusted by shifting values that originally corresponded to a plastic strain of 0.2% to zero, maintaining the value of the yield stress. Adopting the convention of zero plastic strain at the onset of plastic deformation, in line with the principles of continuum mechanics, offers several advantages for ML training. Through this adjustment, model interpretability is improved, and predictions become more intuitive.

Zoom In Zoom Out Reset image size

This shifting process is achieved by applying the following formula:

$$\begin{align}{\boldsymbol{\varepsilon}_p}&= {\boldsymbol{\varepsilon}_{p}}^{MM} - \frac{0.2\% \cdot \boldsymbol{\varepsilon}_{p}^{MM}}{\left\| \boldsymbol{\varepsilon}_{p}^{MM} \right\|} = \left( 1 - \frac{0.002}{\left\| \boldsymbol{\varepsilon}_{p}^{MM} \right\|} \right) \cdot {\boldsymbol{\varepsilon}_{p}}^{MM}\nonumber \\ \left\| \boldsymbol{\varepsilon}_{p}^{MM} \right\| &= \sqrt{\frac{2}{3} \left( \left( {\varepsilon _{p1}}^2 + {\varepsilon _{p2}}^2 + {\varepsilon _{p3}}^2 \right) + 2\left( {\varepsilon _{p4}}^2 + {\varepsilon _{p5}}^2 + {\varepsilon _{p6}}^2 \right) \right)}\end{align} \tag{ 14 }$$

In this equation, ${\boldsymbol{\varepsilon} _{p}}^{{\text{MM}}}$ symbolizes the original plastic strain data from the micromechanical simulations, $\left\| {{\boldsymbol{\varepsilon} _{p}}^{{\text{MM}}}} \right\|$ is the equivalent plastic strain and ${\boldsymbol{\varepsilon} _p}$ represents the adjusted plastic strain data, which aligns with the ML model's assumption of zero plastic strain at the onset of yielding.

Following the preprocessing step and shifting operation, the feature vector takes the form:

$$\begin{align}\left\{ \begin{array}{*{20}{c}} {{x_t} = \left[ {{\sigma _1},{\sigma _2},{\sigma _3},{\sigma _4},{\sigma _5},{\sigma _6},{\text{ }}{\varepsilon _{p1}},{\varepsilon _{p2}},{\varepsilon _{p3}},{\varepsilon _{p4}},{\varepsilon _{p5}},{\varepsilon _{p6}}} \right]} \\ {{y_t} = + 1{\text{ }}} \\ {or} \\ {{y_t} = - 1{\text{ }}} \end{array}\right..\end{align} \tag{ 15 }$$

The corresponding target variable ${y_t}{\text{ }}$ can have two possible values: [−1] or [+1].

The next crucial step in preprocessing is feature normalization. This step is important for many ML algorithms and is even mandatory for SVC to ensure no feature dominates the decision function due to its range. By transforming the features to fall within the range of −1 to 1 it is ensured that all features exert equal influence on the SVC's decision function, which leads to a more balanced and accurate model.

Following the preprocessing step, the prepared training data and corresponding labels can be utilized for model training. Figure 7(a) offers a schematic representation of the processed training data. Here, each yellow point denotes the equivalent stress at yielding under different loading directions. For ease of visualization, only the initial yield surface is shown in a 2-dimensional stress space, as obtained under plane stress conditions. The brown points represent the upscaled plastic data with label +1, while the blue points represent the downscaled elastic data with label −1. It is crucial to highlight that the actual training data are characterized by 12-dimensional feature vectors resulting from the CPFEM simulations in 300 loading directions. The visualization provides only a simplified view of the full dataset. Subsequently, the training process is proceeded with, where the model parameters in the form of support vectors and dual coefficients of the SVC scheme are iteratively adjusted to minimize the discrepancy between the model's predictions and the actual target labels. The process continues until there is either no significant improvement in the model's performance or it reaches an acceptable level of performance. The SVC model employs a kernel function (RBF in this case) to transform the input space, allowing for linear separation of the classes in a higher dimensional space. The visualization in figure 7(b) illustrates how the hyperplane essentially serves as the ML yield function, effectively distinguishing between the elastic and plastic regions within the stress–strain space.

Zoom In Zoom Out Reset image size

Evaluating the performance of a machine learning model is a critical step in assessing its effectiveness and reliability. One commonly used approach is to analyze a confusion matrix, which provides valuable insights into the model's predictive capabilities. The confusion matrix consists of four components: true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN). TP and TN represent the cases correctly identified as positive and negative, respectively. Conversely, FP corresponds to negative cases incorrectly predicted as positive, whereas FN indicates positive cases mistakenly identified as negative. High counts of TP and TN reflect the good predictive power of the model, whereas a high number of FP and FN suggest difficulties in making accurate predictions. From these components, several evaluation metrics can be calculated, including accuracy, precision, recall, and the F1 score. One particularly valuable metric is the Matthews correlation coefficient (MCC), which is well-suited for binary classifications and imbalanced datasets. MCC values range from −1, which signifies a completely inverse prediction, to +1, denoting an ideal prediction. Unlike other metrics, MCC considers every element of the confusion matrix rather than focusing solely on specific components. By penalizing both over-predictions (FP) and under-predictions (FN), it offers a balanced evaluation metric [50]. The MCC score is calculated as follows:

$$\begin{align}{\text{MCC}} = \frac{{{\text{(TP)}}{\text{(TN)}} - {\text{(FP)}}{\text{(FN)}}}}{{\sqrt {\left( {{\text{TP}} + {\text{FP}}} \right)\left( {{\text{TP}} + {\text{FN}}} \right)\left( {{\text{TN}} + {\text{FP}}} \right)\left( {{\text{TN}} + {\text{FN}}} \right)} }} \end{align} \tag{ 16 }$$

It is noted here that measures such as the mean average error are not meaningful to assess the quality of a classification model. Evaluating the performance of an ML model requires the use of well-defined test cases, particularly when employing tools such as the confusion matrix. Ensuring a balanced distribution of positive and negative cases is vital for unbiased performance evaluation. The risk of poorly defined or skewed test cases is that they may result in an unbalanced confusion matrix, thus compromising the integrity and reliability of the derived performance metrics.

To evaluate the robustness and generalization capability of the ML model, it is crucial to utilize unseen data sourced from CPFEM simulations with random loading directions that are not included in the training set. As illustrated in figure 8, these additional test cases are subjected to upscaling and downscaling along their respective hardening paths relative to their associated plastic strain increments. While figure 8 depicts the illustration only for the equivalent flow stress, it is important to emphasize that, as described in section 3.2, this process is applied to the flow stress tensors at each value of plastic strain along the loading path. Through this process, additional test data points are generated that lie within both the elastic and plastic regions of the stress–strain space. These test points follow a Gaussian distribution, clustering densely around the initial flow stress values and gradually decreasing in frequency as the distance from these values increases. The reason for choosing a Gaussian distribution is the fact that points near the boundary are generally perceived as harder to classify, while those farther away are often seen as easier. By this distribution being adopted, a comprehensive and unbiased assessment is ensured, highlighting the adaptability and generalizability of the model when trained across varied scenarios.

Zoom In Zoom Out Reset image size

In this section, the critical preprocessing step necessary for the generation of training data is detailed. Central to this process is the utilization of flow stresses derived from the CPFEM simulations. To create the 12-dimensional feature vectors with corresponding labels, these stresses undergo upscaling and downscaling within a pre-defined data generation range. The discussion explores the determination of the best range for this data generation, ensuring effective coverage of both elastic and plastic regions. In addition, the balance between computational feasibility and predictive performance in deciding the optimal training data size is highlighted. It is important to emphasize that these factors—the data generation range and training data size—are inherently tied to the unique characteristics of the dataset used. Hence, with the introduction of a new dataset, these parameters require a thorough reassessment to ensure optimal performance.

4.1.1. Optimal training data generation range

To find the optimal range for generation of labeled training data, as discussed in section 3.2, our focus was on three critical parameters: ${C_e}$, ${F_e}$, and $N$. Here, $\left[ {{F_e},{\text{ }}{C_e}{\text{ }}} \right]$ is the data generation range, which due to symmetry, is mirrored into both, the elastic and plastic regions of space. The parameter $N$ determines the number of upscaled and downscaled data points in the respective regions. For effective model training, it is important that the training data generation range is optimally selected. Guided by this principle, a low value of 0.65 for ${F_e}$ and a high value of 0.99 for ${C_e}$ were chosen. While this investigation was conducted using $N = 2$, the derived range is applicable for higher $N$ values. A Partial Factorial DoE method was employed to systematically explore specific combinations of values for the parameters ${C_e}$, ${F_e}$. This was done to generate training data within each of these different ranges and to assess their respective impacts of the training result. The selected values for ${C_e}$, ${F_e}$ are summarized in table 4. As explained earlier, the values of ${C_e}{\text{ }}$ and ${F_e}$ represent the distance difference from the actual flow stress considered at position 1 in a 6-dimensional space with suited units for the stress. For example, when ${C_e}{\text{ }}$ is set to 0.85, it implies a spatial distance of 0.15 to the actual flow stress.

Table 4. Range of ${C_e}$, ${F_e}{\text{ }}$used in the Design of Experiment (DoE) for $N = 2$.

Factor	Values
${C_e}$	0.85, 0.9, 0.95, 0.99
${F_e}$	0.65, 0.7, 0.75, 0.8

To simplify the experiment, the data selection was confined to simulations with up to 0.5% plastic strain. ML training was carried out utilizing the training data generated within the range resulting from each combination of the values represented in table 4. The MCC score was computed using designated test cases as detailed in section 3.5, providing a quantitative measure of the trained model's performance under varying conditions, which can be seen in figure 10.

Zoom In Zoom Out Reset image size

Figure 10 shows a heatmap representation of the MCC scores resulting from ML training, where distinct training data sets were generated within different ranges, each originating from a unique combination of ${C_e}{\text{ }}$ and ${F_e}$. Notably, ${C_e}{\text{ }}$ emerges as a more influential parameter. This is evident as, for certain ${C_e}$ values, the MCC scores remain constant or exhibit only minor fluctuations across the entire range of ${F_e}$ levels. Particularly, ${C_e}$ value of 0.9 stands out as optimal, resulting in higher scores than both its immediate neighboring levels. The combination of ${F_e} = 0.7$ and ${C_e} = 0.9$ results in the highest MCC score of 0.8418 and can be chosen as the optimal data generation range.

4.1.2. Optimal training data generation size

As mentioned in section 3.2, during the preprocessing step, for each level of plastic strain, $N$ labeled data points are generated through upscaling the flow stress tensors into the plastic region and another $N$ values through downscaling the flow stress tensors into the elastic region. This process results in a total of $2N$ training data points contributing to constructing the 12-dimensional feature vectors across all plastic strain levels for all the loading directions. The value of $N$ directly dictates the training size, making its optimal selection important. If $N$ is set too high, it can lead to an excessive amount of training data. This not only increases computational costs but can also introduce noise, potentially leading to overfitting and poor performance on unseen data. To determine the optimal value for $N$, ML training was carried out with $N$ values of 2, 4 and 6, using the raw data in the micromechanical stress–strain curves up to 3% plastic strain in the optimal range found for ${F_e} = 0.7$ and ${C_e} = 0.9$. Following training, the MCC scores were assessed using the test cases defined earlier. Along with the MCC scores, the corresponding number of training data and actual support vectors are summarized in table 5.

Table 5. Number of support vectors with respect to that training size in various $N$ and the corresponding ratio.

$N$	MCC Score	Number of Training Data	Number of Support Vectors	Ratio
2	0.8418	157 416	7456	0.047
4	0.8706	314 832	7475	0.024
6	0.8458	472 248	7591	0.016

While there is an improvement in the MCC score when moving from $N = 2$ to $N = 4$, the gain is relatively small. Compared to the significant increase in computational time required for the larger training size, this minor boost in MCC may not be worth the extra computational effort. Additionally, when $N$ is increased to 6, the MCC score drops slightly, which may indicate the start of overfitting as also the number of support vectors is increasing significantly. Support vectors are specific data points within the training set that, during the training procedure, are identified to lie closest to the hypersurface separating elastic and plastic regions in the stress space. As the training size increases, the algorithm has more potential data points available for selection as support vectors. In an ideal scenario, the algorithm must pick the most representative data points as support vectors to define a precise boundary. However, this is not always the case. If the selected support vectors are not representative of the broader dataset's patterns, it typically leads to a decrease in the performance, which can be observed in $N = 6$ case. Given these observations, $N = 2$ emerges as the recommended choice for ML training.

This section sheds light on the results of the ML model trained using the optimal training dataset. Using the optimal training range and setting $N = 2$, resulting in a training size of 157416, the model was trained on data up to 3% plastic strain. A comprehensive hyperparameter optimization process was undertaken with the aim of fine-tuning the model and identifying the optimal combination of hyperparameters. During the optimization phase, a grid search strategy was utilized, facilitating a systematic exploration of a comprehensive range of values for the hyperparameters $C$ and $\gamma$. The values for $C$ were allowed to vary from 1 to 12 while $\gamma$ was explored within the range of 0.1 to 1. To enhance the robustness of this process, 5-fold cross-validation was incorporated, ensuring that diverse subsets of the data were utilized for both training and validation of the model. A comprehensive evaluation was conducted through exhaustive testing of various hyperparameter combinations, allowing for the precise identification of the most effective setting. The optimal combination of hyperparameters, once identified, was employed to train the final model. To assess the model's generalizability and check for overfitting, its performance was evaluated on both the training set and a separate unseen test set within the same range of up to 3% plastic strain, as described in section 3.5. The model's predictive performance is shown in the confusion matrix in figure 11, while the associated performance metrics for both the training and test sets, are given in table 6.

Zoom In Zoom Out Reset image size

Table 6. Performance metrics of the ML model using the optimal hyperparameters on training and test sets.

Set	${\boldsymbol{C}}$	${\boldsymbol{\gamma}}$	Precision	Accuracy	Recall	F1-score	MCC
Train	4	0.5	0.9410	0.9448	0.9492	0.9450	0.8897
Test	4	0.5	0.9216	0.9053	0.8858	0.9034	0.8111

By incorporating plastic strains as additional degrees of freedom, the ML yield function was trained to capture the hardening behavior of the materials. To visualize this behavior, the full 6-dimensional stress space was simplified to the plane of deviatoric stresses in the 3-dimensional space of principal stresses, the so-called $\pi$-plane, which can be conveniently described in polar coordinates comprised of the polar angle $\theta$ and the equivalent stress. Within this plane, the ML yield function is evaluated for different stresses, i.e. points on that plane, and plastic strain tensors, given as additional labels, to assess the yielding behavior under a wide range of conditions. To calculate the proper strain tensor associated to a given level of equivalent plastic strain, the gradient of the yield function is evaluated for each point in the stress space and used to indicate the plastic flow direction. By multiplying this normalized direction with the equivalent plastic strain, any desired strain tensor can be constructed in a realistic manner and provided as argument to the ML yield function. Through the evaluation of the yield function's response for the points on the π-plane, the flow stresses were determined, as shown in figure 12.

Zoom In Zoom Out Reset image size

Figure 12 delineates five distinct contours, each marking the boundary in the stress space where yielding occurs at a specific level of equivalent plastic strain. With each successive contour, the evolution of the yield boundary and the increasing resistance to deformation can be observed. Also, the contours exhibit a smooth, rounded nature, confirming the convexity of the yield surface. Specifically, for any two points selected on a contour, all points on the straight-line segment connecting those two points also lie within the domain bounded by the contour, which is a defining characteristic of convex shapes. This visual representation affirms our expectation of a convex yield function, showcasing a form that is consistent with the physical principles of material behavior when subjected to stress. This means that the ML model captures the essence and primary characteristics of the data without overcomplicating or overfitting. When the training data inherently exhibits a convex yielding behavior, a characteristic expected from the material's physical properties, the SVC model is naturally configured to mirror this convexity. It is seen that the model does not try to introduce convoluted boundaries but aims to replicate the general trend of the data.

In figure 13, the predicted yielding behavior from the ML model is compared with actual scattered data points from the CPFEM simulations, representing two distinct levels of plastic strain: the initial yielding and 2.5% of equivalent plastic strain. For each specified plastic strain level, the flow stresses across all loading directions were collected and then converted into polar coordinates on the π-plane of deviatoric stresses. As shown in figure 13, a remarkable alignment between the predictions of the ML model and the actual empirical data can be observed. Most of the scattered data points are found to nestle close to the contours predicted by the ML function. This strong alignment not only underscores the model's ability to capture the material's intricate hardening behavior but also reaffirms the reliability and robustness of the ML model trained in this study, even in regions where data is relatively scarce.

Zoom In Zoom Out Reset image size

The trained ML model can determine both the initial yielding stress and subsequent flow stresses as plastic yielding progresses. Using this capability, a stress–strain curve for any given loading direction can be reconstructed directly from the trained ML yield function. For demonstration purposes, the micromechanical simulation results from an arbitrarily chosen loading direction used during training were taken. This unit stress tensor, $\sigma = {\text{ }}$ (0.054, −0.091, −0.127, 0.139, 0.55, 0.05), is applied as the boundary condition to the CPFEM model and increased proportionally during the simulation. For each equivalent plastic strain in the dataset, the unique stress multiplier is identified and provided as an argument to the ML yield function. As described above, the corresponding plastic strain tensor has been obtained from the normalized gradient of the ML yield function multiplied by the equivalent plastic strain. The resulting stress–strain curve is depicted in figure 14 and compared to the actual data gathered from the CPFEM simulation. This comparison reveals a significant level of agreement, which underscores the proficiency of the ML yield function in characterizing hardening behavior. It should be noted that the ML model was trained specifically within the range of plastics starting from 0.2%, as discussed above. Therefore, stress–strain data from the elastic branch below this 0.2% threshold were not considered in the plot and would result in a linear-elastic response of the ML yield function.

Zoom In Zoom Out Reset image size