Fast parameter optimization for high-fidelity crystal plasticity simulation using active learning

Material preparation and characterization

A bulk directional solidification sample of Al-Cu alloys with 17.39 at.% Cu was initially prepared in a vacuum induction melting furnace, employing metal raw materials with a purity of 99.99%. A tensile specimen was then extracted from the bulk sample, and its schematic is depicted in Figure 1A, where the long axis is along the withdrawal direction of directional solidification. Macroscopic stretching of the specimen was conducted along the withdrawal direction of directional solidification (x-axis) with an initial strain rate of 10^-4 s^-1 using a TSMT EM6.504 tensile testing machine to obtain its macroscopic stress-strain curve, as shown in Figure 1B. Note that the starting point of the engineering stress-strain curve in Figure 1B does not originate from (0, 0). This is due to the fact that the initial force setting of the force sensor used in the testing procedure typically starts recording data when the force exceeds 54 N. Furthermore, Wei et al. and Ebrahimi et al. have also conducted tensile tests on materials using the same experimental equipment, and their stress-strain curves similarly do not originate from (0, 0)^[34,35]. Their results proved evidence for the reliability of the experimental results obtained in this study. Additionally, it should be noted that all stress-strain curves presented in this study are of the engineering stress-strain curves. Furthermore, an in-situ tensile specimen was fabricated via electro-discharge machining with a gauge length of 8.0 mm, width of 3.2 mm and thickness of 1.0 mm, and its microstructure was observed using a TESCAN MIRA3 field emission scanning electron microscope (SEM), as shown in Figure 1C. The specimen was subsequently ground and finely polished using colloidal silica. A mono-layer of SiO₂ nanoparticles was then deposited on the specimen surface using the drop-casting technique^[36] to ensure a dense and uniform speckle morphology, and Figure 1D displays an example of an emission scanning electron microscope and microscopic-digital image correlation (SEM-DIC) speckle image. The microscopic deformation behavior was observed through in-situ tensile testing, with the specimen being deformed at an initial strain rate of 1 × 10⁴ s^-1 in a Kammrath & Weiss tensile stage. The parameters of the Kammrath & Weiss system can be found at https://www.kammrath-weiss.com/en/tensile-compression-modules/. After testing, the local strain distribution was analyzed using the GOM Correlate software, and details of the local strain distribution will be analyzed in Section “Sensitivity analysis between CP parameters and simulation output”.

Figure 1. (A) Schematics of the uniaxial tensile specimen, where the long axis is the directional solidification direction; (B) The engineering stress-strain curve of the tensile specimen; (C) Microstructure in the directional growth direction of the specimen; (D) An example of SEM-DIC speckle image (C). SEM-DIC: Scanning electron microscope and microscopic-digital image correlation.

CP model

To predict the deformation behavior of Al-Cu eutectic materials, the polycrystalline elasto-viscoplastic CP model proposed by Lebensohn et al. was employed. In this study, we introduced the model briefly, and the details were seen in the literature^[37-40].

For an elastoplastic issue, the total strain for a system $${\epsilon }_{ij}^{\mathrm{t}}(x)$$ can be decomposed into two parts: elastic strain $${\epsilon }_{ij}^{\mathrm{e}}(x)$$ and plastic strain $${\epsilon }_{ij}^{\mathrm{p}}(x)$$. Using Euler implicit time discretization and Hooke’s law, the total strain at point x with t + Δt can be expressed as

where $$C_{ijkl}^{-1}$$ represents the inverse of the elastic tensor C_ijkl. The constitutive relation between the plastic strain rate $${\dot{\epsilon }}_{ij}^{\mathrm{p}}(x)$$ and the stress σ_kl is given by the classic incompressible rate-dependent CP equation:

where $${\dot{\gamma }}^s(x)$$ and $${\dot{\gamma }}_0$$ are the shear rate and reference shear rate, respectively. $$m_{ij}^s(x)$$ is symmetric Schmid tensor with slip system s at point x, n is the stress exponent, $${\tau }_0^s(x)$$ is the critical shear stress (CRSS) of the slip system s, and the symbol “:” represents double contraction of indices.

To solve Equation (2) efficiently, the algorithm of fast Fourier transform (FFT) is employed. In the FFT algorithm, the local deformation response of heterogeneous materials can be expressed by convoluted integrals between the Green function of a linear reference medium and a fluctuating polarization field. Thus, using the stiffness of a linear reference medium and a displacement gradient u_k,l(x), the stress can be expressed as:

here, φ(x) is the polarization field, and the local strain field ε_kl(x) = (u_k,l(x) + u_k,l(x))/2. When the stress equilibrium is reached, σ_ij,j(x) = 0, one can obtain:

As the polarization field described by Equation (5) is a function of an unknown strain field, some iterative procedures are required to solve them. To ensure the effectiveness and stability of solving during iteration, the augmented Lagrangian scheme^[41] proposed by Michel was applied. In this scheme, two auxiliary fields λ and e were introduced, in which $${\lambda }_{ij}^{(i)}$$ and $$e_{ij}^{(i)}$$ represent auxiliary guess stress field and strain fields at the i-th iteration, respectively. The polarization field can be rewritten as:

Thus we obtain:

The augmented Lagrangian scheme is suitable for solving Equation (7), which requires the elimination of a residual (R) for each material point x. The residual R is a function of the stress field σⁱ⁺¹ constitutively related to the strain field εⁱ⁺¹:

Note that, for convenience, in Equation (8), the residual is represented by the contraction symbol for symmetric tensors. For example, σ_ij → σ_k, k = 1, 6; C_ijmn → C_kl, k, l = 1, 6, and so on. The nonlinear Equation (8) is solved using a Newton-Raphson (N-R) method, that is:

Combing Equations (8) and (1), the Jacobian in Equation (9) can be provided as:

The derivative on the right-hand side of Equation (10) is the tangent compliance of Equation (2). Due to the dependence of plastic strain on the stress, the CRSS is a function of stress τ₀^s(ε^p(σ_ij(x)) = τ₀^s(σ_ij(x)). The approximate expression for this tangent compliance is:

Combining Equations (10) and (11), an approximate Jacobian can be given as^[37]:

Despite this approximation, the convergence rate is almost unaffected even in the case of strain hardening, compared to the case of without hardening.

According to the expression of the Jacobian in Equation (11), different hardening laws are obtained. For the strain hardening behavior of Al-Cu eutectics, the CRSS evolution of each slip system is a function of the total accumulated shear strain^[37]:

Where H is a positive scalar related to the strength of the work-hardening behavior. Note that the hardening law employed in this study is phenomenological, and it does not account for dislocation evolution. To encompass more intricate dislocation movements, it is essential to couple the dislocation evolution within the CP model, including dislocation multiplication and annihilation.

Optimization search methodology

Machine learning

Five commonly used machine learning models - linear regression^[42], ridge regression^[43,44], support vector regression (SVR)^[45,46], random forest (RF)^[47,48] and Gaussian Process Regression (GPR)^[49-51] - were used to establish the relationship between CP parameters and mechanical behavior (stress-strain curves). It is important to note that both SVR and GPR models use an isotropic radial-base kernel function (RBF-kernel). The hyper-parameters of the five models were determined by a method of Bootstrap Sampling Cross-Validation (BSCV), and the prediction uncertainty of the models was estimated by 1,000 samples.

The performance of these models was quantified by the root mean square error (RMSE), defined as

(14) $$\text{ RMSE }=\sqrt{\frac{1}{n}\sum _{i=1}^n{({\mu }_i-{\hat{\mu }}_i)}^2}$$

where μ_i and $${\hat{\mu }}_i$$ are the experimental and predicted values at the i-th strain position, respectively. n is the number of discrete points. The smaller the RMSE, the closer the predictions are to experiments.

To accurately elucidate the relationship between CP parameters and a quantitative indicator C_ecp, five candidate machine learning models - Linear, Ridge, SVR, RF and GPR - were trained, and the best model was subsequently selected. Here, C_ecp is used to evaluate the quality of the CP parameter, which lies in characterizing the discrepancy between experimental and simulated stress-strain curves. Further details will be provided in Section “Quantitative evaluation of CP parameter quality”. The hyper-parameter of SVR and RF was determined through exhaustive search and BSCV. The regularization coefficient (C) and the RBF kernel’s length related parameter (γ) were identified as the principal parameters for the SVR model; the number of trees (N) and the maximum tree depth (H) were designated as key parameters for the RF model; the scale parameter(Length_scale) was primarily emphasized within the GPR model; the determination of these hyperparameters is shown in Supplementary Figure 1. The optimal hyperparameters for the SVR, RF and GPR models are listed in Table 1. Figure 2 presents a comparison between the actual and predicted C_ecp by the five models using the initial training set. A significantly larger number of data points deviate from the diagonal line in Figure 2A and B compared to Figure 2C-E. This observation suggests that the Linear and Ridge models fail to capture complex relationships between CP parameters and C_ecp, unlike SVR, RF, and GPR models. Furthermore, nearly all sample points in the GPR model are distributed near the diagonal, underscoring the superior performance of the GPR model relative to other models. Additionally, the RMSE comparison among the five models, depicted in Figure 2F, reveals that the GPR model achieves the lowest RMSE, indicating its highest accuracy in predicting C_ecp. Consequently, the GPR model was chosen to construct the AL framework.

Figure 2. Accuracy evaluation of the trained ML models. (A) Linear regression; (B) Ridge regression; (C) SVR; (D) RF; (E) GPR; (F) RMSE comparison among the models. ML: Machine learning; SVR: support vector regression; RF: random forest; GPR: Gaussian Process Regression; RMSE: root mean square error.

Table 1

Optimal hyperparameters for SVR, RF and GPR models

Models	Hyperparameters
SVR	C = 45, γ = 1.45
RF	N = 95，H = 7
GPR	Length_scale = 1.85

SVR: Support vector regression; RF: random forest; GPR: Gaussian Process Regression.

Search strategy

According to Equations (2) and (13), most CP parameters can be confirmed from Refs^[2,52-55]. Only four parameters remain to be undetermined: $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$, $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$, H_Al and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$. Note that these determined CP parameters come from precise experimental measurements, and can be directly used in CP simulations without the need for calibration. Optimization becomes necessary only for the parameters that are either unavailable from experiments or cannot be experimentally measured, that is $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$, $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$, H_Al and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$. Among them, $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ and $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$ are related to the dislocation movements of α-Al and θ-Al₂Cu, respectively^[37,56,57], while H_Al and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$ are related to the work-hardening behavior of α-Al and θ-Al₂Cu, respectively^[37,57]. Generally, the range of $${\dot{\gamma }}_0$$ is from 10^-2 to 10^-5, with larger values indicating easier dislocation movement. The range of H is from 0 to 1,000 MPa, with larger values indicating stronger work-hardening behavior^[37,57]. If 10^-5 and 10 MPa are the steps of $${\dot{\gamma }}_0$$ and H within their respective ranges, then they have 10³ and 10² possible values, respectively. Thus, the total number of parameter combinations is 10³ × 10³ × 10² × 10² = 10¹⁰. Undoubtedly, searching for the optimal combination in such a large space (denoted as Ω in this study) is a very challenging task. Therefore, we constructed an adaptive framework for CP parameter optimization by borrowing the concept of AL from the field of alloy composition design. The complete process of finding the optimal CP parameters in this study is shown in Figure 3.

Figure 3. Strategy of CP parameter optimization based on AL. CP: Crystal plasticity; AL: active learning.

Part B in Figure 3 illustrates an intelligent framework called AL to guide the iterative refinement of CP parameters. The main process is as follows:

(1) A total of 26 groups of initial parameters were selected in the search space Ω. Using these parameters, 26 simulated stress-strain curves were generated using the CP model, and the corresponding C_ecp values were calculated.

(2) If the maximum C_ecp (denoted as $$C_{\mathrm{ecp}}^{\max }$$) exceeds 0.99, the corresponding inputs are the optimal CP parameters. If not, proceed to the next step.

(3) An appropriate machine learning model f was trained using the initial parameters as inputs and the corresponding 26 C_ecp data as outputs. The model f was then used to predict C_ecp for all parameter combinations in Ω, with the mean and variance of each prediction denoted as μ and σ, respectively.

(4) The efficient global optimization (EGO) utility function was applied to select three parameter combinations by maximizing expected improvement (EI) across the entire Ω. EI is defined as EI = σ[ϕ(z) + zψ(z)], where z = (μ - μ^*)/σ, μ^* is $$C_{\mathrm{ecp}}^{\max }$$, and ϕ(z) and ψ(z) denote the standard normal density and distribution functions, respectively.

(5) The three selected parameter groups were then used as input to the CP model, and their corresponding C_ecp values were calculated.

(6) The initial set was updated with the new recommendations. Steps 2-5 were then repeated until $$C_{\mathrm{ecp}}^{\max }$$ > 0.99, resulting in the identification of the optimal parameter combinations.

It should be noted that the AL strategy above is an adaptive loop that iteratively aligns the CP simulation with experimental results by maximizing C_ecp. Once the optimal parameters are identified, the resulting CP model will have high fidelity. The use of machine learning in this process makes the parameter search intelligent and efficient.

Note that the mechanical properties and microstructure of the Al-Cu eutectic alloy at room temperature were displayed in Figure 1. In the CP modeling, the experimental microstructure serves as the input of the CP model, with comprehensive details provided in Ref^[52]. Additionally, the initial configurations and material parameters of the CP simulation for the Al-Cu eutectic materials are from Ref^[52].

Quantitative evaluation of CP parameter quality

The key to evaluating the quality of the CP parameter lies in characterizing the discrepancy between experimental and simulated stress-strain curves. To quantify this discrepancy, we discretized both the experimental and simulation curves into 17 data points, ensuring that the horizontal coordinates of each set of points were consistent. Then, the determination coefficient R² was used to assess the consistency between discrete experimental and CP simulation data. Notably, to distinguish it from the accuracy evaluation of machine learning illustrated in Figure 2, we redefine it as C_ecp, that is:

(15) $$C_{\mathrm{ecp}}=R^2=1-\frac{\sum _{i=1}^n{({\mu }_i-{\hat{\mu }}_i)}^2}{\sum _{i=1}^n{({\mu }_i-\overline{\mu })}^2}$$

where $${\overline{\mu }}_i$$ represents the average of n discrete points. The closer C_ecp is to 1, the closer the simulations are to the experiments.

To validate the effectiveness of C_ecp in assessing CP parameter quality, Figure 4 shows the degree of alignment between two examples of simulated stress-strain curves (simulation 1 and simulation 2) and the experimental curve. The CP parameters of the two examples are listed in Table 2. Observations from Figure 4 revealed that, in the nearly linear portion of the stress-strain curve, the data of Simulations 1 and 2 are in good agreement with the experimental data, demonstrating that the elastic parameters of the simulations align with those of experiments. This concurrence is attributed to the fact that the elastic parameters utilized in this study come from experimental measurements. Upon departing from the initial linearity, the data points from Simulation 2 (C_ecp = 0.65) significantly deviate more from the experimental data than those from Simulation 1 (C_ecp = 0.86). This discrepancy implies that the dislocation movements (reflected by parameters $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ and $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$) and the work-hardening behavior (reflected by parameters $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ and $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$) in Simulation 2 are more deviated from the experiment than that in Simulation 1. Additionally, the curve of simulation 1 is closer to the experimental curve than that of simulation 2. The C_ecp value between simulation 1 and experiment is 0.86, which is higher than that for simulation 2 (0.65), as shown in Table 2. This indicates that a higher C_ecp value corresponds to a smaller disparity between simulation and experiment. This further demonstrates that C_ecp is indeed a valid metric for quantifying the quality of CP parameters. In this study, 0.99 is selected as a threshold of C_ecp.

Figure 4. Example of quantitative evaluation of CP parameters. CP: Crystal plasticity.

Table 2

Qualitative and quantitative assessment of the degree of agreement between simulation and experiment under two sets of parameters

	$${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$	$${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$	HAl (MPa)	$$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$ (MPa)	The degree of agreement	Cecp
Simulation 1	0.01	0.00001	10	100	Proximity	0.86
Simulation 2	0.007	0.0035	400	900	Deviation	0.65

Iterative optimization of CP model parameter

At the early stage, a trial-and-error approach was performed to establish an initial training set, as detailed in the Supplementary Materials. The maximum C_ecp in the initial training set is 0.95, indicating that the simulated stress-strain curve deviates far from the experimental curve. Following the AL framework shown in Figure 2, we developed a high-precision GPR model, which was utilized to iteratively improve the CP parameters. It is important to note that when using a GPR-based search strategy to recommend new parameter combinations, the constraint must be met that the parameter $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$ must be less than the parameter $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$. This is because, during the deformation of Al-Cu eutectic alloys at room temperatures, the activation of dislocations in the θ-Al₂Cu phase is more difficult than that in the α-Al phase.

As illustrated in Figure 5, C_ecp surpassed 0.99 in just four iterations, and a total of six iterations were executed. Detailed recommendations are provided in Table 3. In the first iteration, there is a significant difference among the three points, and the maximum C_ecp is only 0.86, which is below the threshold. In the second iteration, the C_ecp values noticeably increase compared to those from the first iteration, and the data range significantly reduces. This indicates that the GPR samples are centrally located at a specific position of data sparsity in Ω. Notedly, C_ecp continues to rise until it reaches 0.99 in the fourth iteration. Interestingly, in this iteration, the three sets of CP parameters are very close to each other, suggesting that the global optimal may be successfully identified. As there is no more suitable data in Ω to further improve C_ecp, its values start to decrease from the fifth iteration onward. By the sixth iteration, the data points began to disperse again, indicating that the exploration in Ω guided by GPR has resumed to improve the data distribution of the training set. Furthermore, the C_ecp values of the recommended CP parameters in the 5th and 6th iterations are consistently lower than the C_ecp value obtained after the 4th iteration, indicating that the optimal CP parameters have been identified. Notedly, C_ecp values in the 5th and 6th iterations have been over the preset threshold of 0.99 but the iteration continues, which is just to make it more intuitive that the CP parameters in the 4th iteration are optimal. Additionally, as shown in Table 3, among the set of optimal CP parameters that correspond to the maximum C_ecp value, the parameter H_Al is approximately 100 MPa and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$ is around 300 MPa. This consistency with the sensitivity analysis presented in Section “Comparative analysis between the AL strategy and the typical traditional methods” further confirms the discovery of the optimal CP parameters. Since the optimal has been discovered, the iterative search is terminated.

Figure 5. C_ecp for the recommendations in a total of six iterations. Iteration 0 corresponds to the initial training data.

For an AL process that maximizes EI, the initial sparser training data leads to greater prediction uncertainty of GPR on the new sample points and increases the randomness of data distribution^[58], as shown in the first iteration in Figure 5. The prediction accuracy on each recommendation is depicted in Figure 6A. Evidently, GPR overestimates the C_ecp for the three new sets of CP parameters, resulting in significant prediction errors and uncertainties. From the mathematical expression of the utility function EI that is used to recommend the next experimental points, it is easy to understand that given a large uncertainty σ, EI → σ. That is, in the limit of high uncertainty, EI is maximized by choosing the next best material with the largest uncertainty (“exploration”). It represents a region in our search space, where the regressor is most uncertain. Thus, for regions not covered by the training data, GPR will also centrally sample in the early stages. Consequently, the values of C_ecp at the second iteration are almost consistent, but they deviate significantly from the predictions, as indicated by the orange points in Figure 6A. Once the training set’s sparse area is filled, prediction uncertainty decreases, model accuracy increases, and the search moves closer to the optimal solution^[58]. On this occasion, EI → μ - μ^*, where μ is the predicted average by the machine learning (ML) model, μ^* is the current maximum value in the training set. This means that in the limit of no uncertainty or if uncertainties are almost the same for all points, EI is maximized by choosing the largest μ (i.e., “exploitation”). Therefore, the next best material for measurement is chosen by faithfully trusting the predictions from the regressor. This conclusion is supported by the trends of the green and red dots (at the fourth and fifth iterations) in Figure 6A. After identifying the global optimization, GPR explores other sparse regions in Ω, leading to a more dispersed distribution of new sampling points, yet with accurate predictions, as shown by the purple and brown spots in Figure 6A. Additionally, Figure 6B shows that EI rapidly decreases from the first to the fifth iteration. Once the maximum C_ecp has been identified, the probability of new recommendations improving upon the best training data is reduced to 0, causing EI to drop to 0.

Figure 6. Prediction accuracy and EI variation during the optimization process of CP parameters. (A) Accuracy between actual and predicted C_ecp; (B) Variation of EI with the number of iterations. EI: Expected improvement; CP: crystal plasticity.

Mechanical behavior analysis of Al-Cu eutectic materials

To validate the optimal CP parameters, we next input the parameter into the CP model, obtained a simulated stress-strain curve, and compared it with the experimental curve. It is crucial to emphasize that, in this study, the input microstructure for the CP model is firmly based on direct experimental observations^[52]. For quantitative comparisons between the experimental data and the CP simulations, the horizontal axes of the experimental and simulated stress-strain curves have been aligned. For comparison, other curves for the best parameter combinations at each AL iteration are also included, as shown in Figure 7. Note that the CP simulation of a complex stress status with the xy shear deformation has been provided in Supplementary Figure 2 to ensure the validation of the CP simulation. From Figure 7, it is evident that the simulated stress-strain curve for the second iteration (orange curve) is closer to the experiment (black curve) than that for the first iteration (blue curve). This is also reflected in the comparison of the C_ecp list in Table 3. Similar patterns can be observed in other curves, demonstrating that C_ecp is a suitable quantitative index for evaluating the reliability of CP simulation. Notably, the red curve in Figure 7 nearly coincides with the black dots, indicating that the recommendations from the fourth iteration may be optimal. In addition, the deviation of the purple (fifth iteration) and brown (fifth and sixth iterations) curves in Figure 7 from the black dots indicates that the recommended parameters begin to deviate from the optimal solution after the fourth iteration. This further demonstrates the recommendations from the fourth iteration are indeed optimal, resulting in a high-fidelity simulation that closely matches reality. Note that the elastic part of CP simulated stress-strain curves with different iterations is consistent due to the fact that the elastic parameters in this study are sourced from literature^[55,59,60].

Figure 7. Comparison of simulated and experimental stress-strain curves. The simulations are performed with the best parameter combinations at each AL iteration. AL: Active learning.

Figure 8 shows a comparison between the simulated local strain distributions and the results from in-situ DIC experiments. The CP simulation was performed using the optimal parameters. As shown in Figure 8A(i)-(iii), with rising total strain ε_t, experimental strain concentration at microstructural defects increases, marked by the black dashed lines. Similarly, in the CP simulations, as shown in Figure 8B(i)-(iii), an analogous trend in strain concentration at microstructural defects with increasing total strain is observed, also marked by the black dashed lines, demonstrating the consistency of the experiments and simulations in microscopic deformation. This further illustrates the reliability of the optimal parameters calibrated using our proposed search strategy. In addition, when the total strain is 0.001 or 0.0025, the minimum strain is not equal to 0 in CP simulation, as shown in Figure 8B(ii) and (iii), which is different from the experimental observation shown in Figure 8A(i) and (iii). The discrepancy arises from the fact that the CP model does not account for dislocation reactions between different slip systems, such as dislocation annihilation and multiplication. In the actual room-temperature tensile process of the Al-Cu eutectic material, there exist dislocation reactions (dislocation multiplication or annihilation) in the α-Al phase, and more significantly, these mobile dislocations can cross the two-phase interface and enter the Al₂Cu phase, further affecting the local strain within the Al₂Cu phase^[52,61-63]. However, since our CP model ignored the dislocation multiplication or annihilation, there is a discrepancy between the simulation and experimental results as depicted in Figure 8. It is worth noting that although our CP model primarily investigates the situation where the Al₂Cu phase is without plastic deformation, it remains applicable to cases where considerable plastic deformation is present in the Al₂Cu phase. For further details, refer to Supplementary Figure 3.

Figure 8. Comparison between the simulated and experimental local strain distributions. The CP simulation was completed with the optimal parameters. CP: Crystal plasticity.

To further explain microscopic deformation behavior, we examined the stress and plastic strain distributions at different total strains, as shown in Figure 9. In Figure 9A(i)-(iii), the red phase and blue phase represent α-Al phase and θ-Al₂Cu phase, respectively. Figure 9B(i)-(iii) shows that the minimum stress ranges from 4 to 8 MPa, while the maximum stress increases from 8.5 to 46 MPa as the total strain ε_t increases from 0.0005 to 0.0025. In Figure 9C(i)-(iii), the maximum plastic strain gradually increases, whereas the minimum plastic phase remains at 0. The minimum stress and the maximum plastic strain are observed in the α-Al phase, while the maximum stress and minimum plastic are found in the θ-Al₂Cu phase. This indicates that plastic deformation occurs in the α-Al phase even under small stresses, whereas the θ-Al₂Cu phase is subjected to large stresses without plastic deformation.

Figure 9. Phase field, stress, and plastic strain distributions with different total strains in CP simulations. (A) (i-iii) Phase field; (B) (i-iii) Stress distributions; (C) (i-iii) Plastic strain distributions. CP: Crystal plasticity.

To further support this conclusion, we analyzed linear stress and plastic strain distributions along the x-axis [labeled by the black dotted line in Figure 9B(i)], and the results are shown in Figure 10. It was found that as the total strain ε_t increases, the stresses in both the α-Al and θ-Al₂Cu phases increase, with the stress in the θ-Al₂Cu phase being significantly higher than that in the α-Al phase. As the total strain increases, the increased rate of the stress difference between σ_max and σ_min in the α-Al phase is significantly slower than that in the θ-Al₂Cu phase, as shown in Table 4. This indicates that the stress variation in the θ-Al₂Cu phase is much greater than that in the α-Al phase. Moreover, only the plastic strain in the α-Al phase shows a gradual increase, while plastic deformation in the θ-Al₂Cu phase remains at 0. For dual-phase materials, the deformation behavior is not only related to the stress and strain distributions within each phase but also significantly influenced by the inhomogeneous deformation near the interfaces between the two phases. In this study, we quantified the non-uniform deformation near the α-Al and θ-Al₂Cu phase interface under different total strains. To better investigate the inhomogeneous deformation near the phase interface, the maximum plastic strain difference of the α-Al and θ-Al₂Cu phases near the interface with different total strains was counted. Since the plastic deformation of the θ-Al₂Cu phase is zero, the inhomogeneous deformation of the α-Al and θ-Al₂Cu phase dual-phase material is primarily affected by the α-Al phase. The statistical results show that under total strains of 0.0005, 0.001, and 0.0025, the maximum plastic strain in the α-Al phase is 0.000018, 0.000375, and 0.00317, respectively. This indicates that as the total strain grows, the inhomogeneous deformation of the Al-Cu eutectic alloy significantly increases. These findings are not only consistent with the DIC experiment of Figure 8 but also align with the previous literature reports^[52,62-64], further confirming the reliability of our search strategy in optimizing CP parameters.

Figure 10. Line stress and plastic strain distributions with different total strains. (A) Line stress; (B) Plastic strain distributions. The location of the lines is labeled by the black dotted line in Figure 9B(i).

In summary, the AL-based search strategy successfully optimized the CP parameters for the Al-Cu eutectic alloy at room temperature. Although the search strategy in this study was only applied to the parameter optimization of a quite simple CP model with only four parameters, it can be anticipated that the strategy proposed in this study is also suitable for more complex systems, such as torsional deformation^[65,66]. In this case, similar treatment for the parameters and the quantification of the simulation reliability should be performed, and more iterations may be required due to more undetermined parameters than those in our study, but even so, it is faster than the trial-and-error selection. However, we also acknowledge that for the parameter optimization of non-uniform computations^[67,68], our strategy may need to be further improved.

Sensitivity analysis between CP parameters and simulation output

To quantify the impact of the CP input parameters ($${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$, $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$, H_Al and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$) on the CP simulation output C_ecp, we assumed four cases by a control variable method. The predicted results by GPR are given in Figure 11. Note that four cases in this study were determined randomly yet all meet the restriction that $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ is greater than $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$ due to the slip resistance of α-Al phase being significantly lower than that of θ-Al₂Cu phase at room temperature^[2,53].

Figure 11. Sensitivity analysis between CP parameters and C_ecp. (A) $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$vs. C_ecp with a constant of $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$, H_Al and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$; (B) $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$vs. C_ecp with a constant of $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$, H_Al and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$; (C) H_Alvs. C_ecp with a constant of $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$, $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$ and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$; (D) $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$vs. C_ecp with a constant of $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$, $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$ and H_Al. CP: Crystal plasticity.

Figure 11A displays the variation of C_ecp with respect to $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ while $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$ (0.00001), H_Al (500 MPa) and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$ (500 MPa) are held constant. It is evident that as $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ increases, C_ecp decreases linearly. This is due to the fact that the dislocation movement in the α-Al phase is increasingly arduous with the increasing parameter $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$, resulting in elevating the yield strength of the α-Al phase. Considering that the yield strength of the α-Al phase is inferior to that of the Al-Cu eutectics, a reduction in $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ brings the CP simulated outcomes closer to the experimental data, that is a larger C_ecp. Similarly, as the parameter $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$ diminishes, the dislocation movement in the θ-Al₂Cu phase is increasingly hindered, rendering CP simulation outcomes more congruent with the experimental findings (a greater C_ecp), as shown in Figure 11B. Herein, $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ (0.01), H_Al (500 MPa) and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$ (500 MPa) are maintained as fixed values.

Figure 11C displays the variation of C_ecp with respect to H_Al while $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ (0.01), $${\dot{\gamma }}_{0\mathrm{\_}{\mathrm{Al}}_2\mathrm{Cu}}$$ (0.001) and $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$ (500 MPa) are held constant. As H_Al increases, an initial rise is followed by a decline in C_ecp, and the peak value appears when $${\dot{\gamma }}_{0\mathrm{\_}\mathrm{Al}}$$ is near 100 MPa. H_Al is related to the hardening behavior of phase α-Al, and the hardening behavior of the material remains essentially unchanged under a given deformation condition (room temperature). Therefore, with the variation of the H_Al parameter, an optimal point of H_Al exists at which C_ecp attains its maximum value. Similarly, regarding the $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$ parameter, its correlation with the hardening behavior of phase θ-Al₂Cu is noted, and a comparable trend between $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$ and C_ecp is observed, as shown in Figure 11D. Notably, the peak C_ecp value is observed when $$H_{{\mathrm{Al}}_2\mathrm{Cu}}$$ is around 300 MPa.

From the analysis above, it can be concluded that the AL-based CP parameter optimization strategy has effectively quantified the influence of CP input parameters on simulation reliability for Al-Cu eutectic alloys. This strategy can be easily extended to the optimization of other CP parameters, in which the parameters required to be optimized must lie within a certain range, and a quantitative indicator should be employed to evaluate the deviation between CP simulation results and the experiments. For example, as reported in Ref^[2], the parameters $$g_0^s$$, h, $${\dot{\gamma }}_0$$ and m can also be determined by our proposed strategy. Employing these parameters as inputs of the strategy and the discrepancies between simulations and experiments as the output, one can rapidly determine the most optimal group of these parameters in several iterations.

Comparative analysis between the AL strategy and the typical traditional methods

The traditional methods used for optimizing CP parameters mainly include trial-and-error and inverse analysis. The former is the most common and the simplest method. It involves randomly selecting an initial set of parameters to construct a CP model and evaluating the discrepancy between CP simulations and experimental results. This process is repeated until either an acceptable discrepancy is achieved or computing resources are excessively consumed^{[2,6,8,11-13]}. Notably, the discrepancy is often qualitatively judged by visual observation of stress-strain curves or other patterns. Additionally, the arbitrary nature of determining when to stop CP computation can lead to excessive cost consumption^[69,70]. Moreover, the results of the trial-and-error method depend on the operator’s experience and intuition. When a large number of parameters need to be determined, there is a significant risk of not quickly finding the optimal solution. Reviewing the AL strategy proposed in Figure 2, we observe a distinct advantage: it quantifies the discrepancy between simulated and experimental stress-strain curves by a simple index, C_ecp, whose reliability has been demonstrated in Figures 3 and 7. Notedly, C_ecp is mathematically identical to R². The renaming to C_ecp is to distinguish it from R² for the performance evaluation of the machine learning algorithms used in this study. The usage of C_ecp enhances the comprehensibility of our methodology and simultaneously ensures the ease of data handling. Moreover, C_ecp provides a quantitative bridge that guides the simulation process from experimental data. Furthermore, by using our strategy, as shown in Figures 5 and 6, and Table 3, intelligent searching of high-dimensional parameters has been realized, which significantly reduces the computational overhead in terms of time and cost, requiring only six iterations.

The inverse analysis method begins by proposing a set of initial parameters, and then evaluates the deviation between CP simulation and experimental results within the trust region by a trust-region policy optimization algorithm^[71,72]. If the deviation is large, the algorithm adjusts the trust region, repeating the process until the deviation is minimized. The inverse analysis allows for a direct connection between parameters and simulated results (e.g., mechanical behavior of materials), but it can lead to unnecessary resource wastage^[73]. Additionally, the mathematical modeling of the algorithm in the trust region is highly complex, and the results significantly depend on the rationality of the initial guessed parameters^[73]. Generally, optimization algorithms are also needed to ensure the uniqueness of the solution^[73]. Unlike the inverse analysis, our proposed strategy achieves optimization by maximizing the EI of candidate parameter combinations. This approach balances “exploration” (improving data distribution) and “exploitation” (approaching the optimal) within only a few iterations^{[20,23,25,26]}. The search results are unaffected by the initial dataset as the model input is updated from each iteration. Although we use a non-parametric model, GPR, to ensure accuracy and efficiency, the choice of the surrogate models in our strategy is flexible, including other ML algorithms and genetic algorithms.

We acknowledge that this search strategy has limitations and may only be suitable for some specific scenarios, such as when the material’s deformation behavior is adequately represented by stress-strain curves or other one-dimensional discrete data. Note that the CP parameters in this study were optimized by a single stress-strain curve from an experimental measurement. However, for higher accuracy, stress-strain curves should be based on multiple replicate experiments. So, to better optimize the CP parameters, it is necessary to optimize them using data from multiple replicate experiments in future work.