A virtual laboratory based on full-field crystal plasticity simulation to characterize the multiscale mechanical properties of AHSS

In this work, we proposed a virtual laboratory based on full-field crystal plasticity (CP) simulation to track plastic anisotropy and to calibrate yield functions for multiphase metals. The virtual laboratory, minimally, only requires easily accessible EBSD data for constructing the highly-resolved microstructural representative volume element and macroscopic flow stress data for identifying the micromechanical parameters of constituent phases. An inverse simulation method based on a global optimization scheme was developed to identify the CP parameters, and a nonlinear least-squares method was employed to calibrate yield functions. Mechanical tests of advanced high strength steel sheet under various loading conditions were conducted to validate the virtual laboratory. Three well-known yield functions, the quadratic Hill48 and non-quadratic Yld91 and Yld2004-18p yield functions, were selected as the validation benchmarks. All the studied functions, calibrated by numerous stress points of arbitrary loading conditions, successfully captured both the deformation and strength anisotropies. The full-field CP modeling correlated well the microscopic deformation mechanism and plastic heterogeneity with the macromechanical behavior of the sheet. The proposed virtual laboratory, which is readily extended with physically based CP models, could be a versatile tool to explore and predict the mechanical property and plastic anisotropy of advanced multiphase materials.


Theory
Finite strain CP theory. For completeness, the adopted phenomenological CP theory is briefly presented; more details can be found in Zhang et al. 28 . The CP model was implemented in a finite strain framework based upon the classical multiplication decomposition of deformation gradient as follows, where F p describes the plastic slip on slip planes and F e describes the elastic stretching and the rotation of lattice. The plastic velocity gradient is expressed as, where γ α is the slip rate of the α-th slip system, and S α is the Schmid tensor. m α and n α denote the slip direction and the normal of the α-th slip system. The slip rate γ α is described by a phenomenological plastic flow law as follows, where γ 0 is the reference shear strain rate, τ α the resolved shear stress acting on the slip system, and T e the second-order Piola-Kirchhoff stress tensor. m represents the strain rate sensitivity and g α is the slip resistance that increases asymptotically from the initial value g 0 to the saturation value g s . To represent the work hardening behavior of crystalline materials, the evolution of g α is formulated by a rate-form hardening model as follows 29 , www.nature.com/scientificreports/ where h 0 , g ∞ , and a are material parameters representing the reference self-hardening coefficient, the saturation value of the slip resistance, and the hardening exponent, respectively. Besides, the latent hardening parameter q is routinely assumed to be 1.0 for coplanar slip systems and 1.4 for non-coplanar slip systems 30 .
With the assumption of hypo-elasticity of most metals, T e can be calculated from the elastic Green strain tensor E e as where C is the elasticity tensor; for materials with a cubic crystal structure, C is specified by three parameters, i.e., C 11 , C 12 , and C 44 .
The above CP-based constitutive model in conjunction with an open-source spectral solver DAMASK 31 enables the VL to carry out a large number of virtual experiments with low computational cost 23 . Modeling setup for full-field CP simulations. The geometrical model used for the full-field CP simulations and virtual tests, i.e., the multiphase microstructural representative volume element (RVE), was constructed via the open-source code DREAM.3D 32 with the processed EBSD data as input. The cubic RVE, as shown in Fig. 1, contains the essential microstructural characteristics of the DP980 sheet, including the volume fraction, grain size, and ODF of both ferrite and martensite phases. The RVE contains 100 elements (the Fourier grids with grid size of 0.5 μm) in each direction and has 872 orientations ("grains") in the ferrite phase and 4193 orientations ("islands") in the martensite, respectively. Periodic boundary conditions were applied to the RVE to perform the simulations. For the simulations of uniaxial tension in different directions from RD to TD, for the convenience of calculating r-values, the method of rotating Euler angles of grain orientations instead of alternating the load directions is adopted; more details can be found in Ref. 22 . Following the Bunge convention of Euler angles, the orientation after an in-plane rotation from RD to TD at angle θ is {ϕ 1 − θ , �, ϕ 2 }. In this way, the same RVE and PBCs were used for the full-field CP simulations, and a large number of random stress states were simulated for calibrating the phenomenological yield functions.
Parameter identification for the CP model. The dual-phase steel has a composite microstructure of relatively strong martensite islands dispersed in the soft and ductile ferrite matrix. Both the ferrite and the martensite phases have a BCC crystal structure (herein, we ignore that slight distortion of the martensite's structure from the BCC structure) with the most easily activated slip system families 1 10 �111� and 2 11 �111� at room temperature; each slip system family has an individual set of parameters. Thus, there are a lot of material parameters to be identified. To identify the mechanical properties of individual phases in AHSS, Hu et al. 33,34 employed the state-of-the-art high-energy X-ray diffraction technique and in-situ tensile tests to measure the lattice strain of each phase; they calibrated the CP parameters by correlating the predicted lattice strains with the experiments. This approach is attractive and persuasive but not easily accessible. In this work, an inverse simulation method evolving global optimization was developed to identify the CP constitutive parameters of the ferrite and martensite phases. First, the elastic constants C 11 , C 12 , and C 44 , the reference slip rates γ 0 , and the rate sensitivity coefficients m were assumed identical for the two slip systems of each phase, and their values were routinely documented 35 . It is noted that the elastic properties of the constituent phases in multiphase steels are not fully understood. Souissi and Numakura 36 gave a comprehensive study on the elastic properties of two types of martensite. The four parameters, g 0 , g ∞ , h 0 , and a , which describe the work hardening behavior of slip systems, were identified by an in-house inverse simulation procedure based on the uniaxial tensile data in RD and TD. According to the research of Tasan et al. 35 and Hamada et al. 37 , the ratio of g 0 in the 1 10 �111� and  www.nature.com/scientificreports/ 2 11 �111� slip systems is approximately 1.0 ~ 1.1. Thus, a constraint was enforced to make the initial and saturation slip resistances g 0 and g ∞ of system 1 10 �111� smaller than those of 2 11 �111�.
The inverse simulation method for parameter identification consists of a global optimization platform and full-field CP simulations. The particle swarm optimization (PSO) toolbox of MATLAB® (version R2018b) was employed to identify the global optimum material parameters for the ferrite and martensite phases. In the beginning, reasonable bounds of the parameters were provided for the inverse simulation platform. The platform generated different initial guesses for the parameters, launched numerous full-field CP simulations of uniaxial tension parallelly, and finally extracted the predicted true stress-strain data. The obtained true stress-strain data were compared with the experiment ones to evaluate the root-mean-square deviations (RMSD) R c of the objective function in the current generation. Then the PSO toolbox introduced a new generation with optimized parameter sets. The fitting process will be iterated until the minimum R c below the prescribed limit. Figure 2 shows the flow chart of the inverse simulation procedure for identifying the CP constitutive parameters. The objective function of the PSO optimization procedure reads, where x = [x 1 , x 2 , . . . , x n ] are the parameters to be identified. σ sim i and σ exp i are the simulated and experimental stresses at the same strain, respectively, and i denotes the corresponding data in different directions, i.e., RD and TD. w i is the weighting factor, which is set to 1.0 in this work. Table 1 lists all the CP constitutive parameters used in this work. Figure 3 shows the experimental and the CP simulated flow stress curves of uniaxial tension in both RD and TD. For g 0 , the values of two slip systems of 1 10 �111� and 2 11 �111� in the martensite are 4 ~ 5 times higher than those in the ferrite; for g ∞ , the values of 1 10 �111� and 2 11 �111� in the martensite are 3 ~ 4 and 1.5 ~ 2.5 times higher than those in the ferrite, respectively.
Virtual laboratory. The VL is executed automatically with a Python script. It launches enough full-field CP simulations of the constructed RVE under arbitrary loading conditions generated by a random number generator, post-processes the simulation results, extracts the yield stress points (stress tensors), and then identifies the yield functions. As the illustrative benchmarks, three well-known yield functions, i.e., the quadratic Hill48 5 and Figure 2. The flow chart of the in-house inverse simulation method for identifying the CP model parameters. www.nature.com/scientificreports/ non-quadratic Yld91 9 and Yld2004-18p 12 yield functions, were calibrated to predict the anisotropy of the DP980 sheet.
The nonlinear least-squares method (NLSM) was adopted to identify the coefficients of the studied yield functions. The NLSM problem is solved by using a bounded Levenberg-Marquardt optimization algorithm 23 . With the analytical Jacobian matrix, this algorithm provides enough robustness to find the solution of strongnonlinear issues, even though the initial guess starts far off the final minimum. The objective function O(β) is defined as follows.
φ σ i , β j denotes the effective stress calculated by the calibrated yield functions, and σ i (i = 1, 2, . . . , N) are the set of stress tensors of the corresponded virtual tests at the same plastic work per unit volume (shorten as the specific plastic work in the following content), and β j = β 1 , β 2 , . . . , β M are the parameters of the yield functions to be calibrated. The optimal parameters β j minimize the objective function O(β) . All the stress tensors were assigned with the same unit weight in this study; however, different weights can be endowed to realize the strong dependence of certain stress states. To stabilize and to accelerate the solving process of the objective function  www.nature.com/scientificreports/ O(β) , analytical Jacobian matrices were derived for the yield functions; the form of the Jacobian matrices reads as follows, with N and M denote the number of stress points and the number of material parameters, respectively. For complex yield functions, e.g., the Yld2004-18p, the chain rule was employed to get the Jacobian matrices. More details can be found in the Supplementary. The fitting quality is measured via RMSD ( R y ) of the residual as follows, Apart from the VL, which was used to generate the random stress points, full-field CP simulations of uniaxial tension at a 7.5 o interval from RD to TD were also executed to predict the in-plane anisotropy of the DP980 sheet, and the same RVE was used. Finally, r−values (the Lankford coefficients) and normalized yield stresses were extracted from the results of the CP simulations.

Results and discussion
Microstructure. The studied material is a commercial cold-rolled and annealing DP980 steel sheet with a thickness of 1.2 mm provided by Baosteel. The chemical composition of the steel is summarized in Table S1 in the Supplementary. Figure 4 displays the microstructure of the material characterized on the RD-TD plane in the center of thickness. Figure 4a is the orientation imaging map (OIM) obtained from EBSD data, and the OIM is colored with the RD inverse pole figure (IPF). Based on the band contrast value, as shown in Fig. 4b, the original microstructure was separated into the individual ferrite and martensite phases 17 . Figure 4c presents the orientation distribution function (ODF) maps reconstructed from the EBSD data of two individual phases. Both phases show the typical texture of BCC metals after cold-rolling and annealing operations, i.e., the γ-fiber consisting of {111}�uvw� orientations and the α-fiber consisting of {hkl}�110� . However, the dominant texture components in the ferrite phase are different from those in the martensite phase. For the ferrite, the intensity of α-fiber is quite smaller than that of γ-fiber, and γ-fiber indeed concentrates at the texture component of {111}�121� ; for the martensite, on the contrary, the intensity of α-fiber is comparable with γ-fiber, and α-fiber concentrates at {001}�110� and the γ-fiber at {111}�121� . It suggests sufficient annealing of the ferrite phase, as the annealing process of cold-rolled BCC materials strengthens γ-fiber and reduces the intensity of α-fiber 22 .
The flow stress directionality and yield loci of the DP980 sheet. Figure 5 presents the flow stress curves of the DP980 sheet subjected to different loading conditions. Figure 5a shows the curves of uniaxial tension in different directions, and Fig. 5b the curves of pure shear and biaxial tensile tests. The results exhibit the typical characteristics of cold-rolled dual-phase steels, i.e., no obvious yield point and no stress plateau but continuous hardening. The von Mises equivalent stress of pure shear is higher than the flow stresses of uniaxial tension in RD and biaxial tension with different stress ratios, and the biaxial tensile flow stresses are a little smaller than that uniaxial tensile one in the early deformation stage. With the increase of deformation, the equibiaxial flow stress is higher than the uniaxial tensile stress. Figure 6 presents the 2D yield loci of the calibrated yield functions in the σ 11 − σ 22 space, 60 stress points (open circles) randomly generated by the VL, and experimental results (solid squares) of uniaxial tension in RD and TD, pure shear, and biaxial tension with different stress ratios. For the studied DP980 sheet, the equivalent equibiaxial tensile stress is almost identical with the uniaxial tensile stress in RD but slightly higher than that in TD. The VL simulated and experimental stress points were selected from the deformation stages with the specific plastic work of 3 MPa (Fig. 6a) and 8 MPa (Fig. 6b). All the stress points were normalized by the corresponding uniaxial tensile stress in RD. To verify the VL's capacity in predicting the plastic anisotropy of multiphase materials, only the simulated 60 stress points and no r-values were used to calibrate the yield functions. Table 2 lists the identified parameters of these functions. Different from the recommendation of Hosford 7 , i.e., m = 8 for materials with FCC crystal structure and m = 6 for BCC structure, the intentionally calibrated homogeneous exponent m varies in the range of 5 ~ 8.
As shown in Fig. 6, the yield loci predicted by the VL (the randomly distributed stress points represented by the open circles) agree well with the experimental one (enclosed by the solid squares), i.e., the VL successfully reproduced the yield loci of the DP980 sheet. All the calibrated yield functions finely outline the yield loci of the sheet in comparison with both the experimental and the VL results. In particular, the yield loci predicted by the Yld91 and Yld2004-18p functions are quite coincident. In contrast, those predicted by the Hill48 yield function show a large deviation from others in the regions close to equibiaxial tension and plane strain. This is because the Hill48 yield function is insufficient in capturing the large curvature change of yield locus in the region close to equibiaxial tension. Figure 7 shows a direct comparison of yield loci of all the considered yield functions at equal normalized shear stress σ 12 /σ . The maximum σ 12 /σ of the Yld2004-18p, Yld91, and Hill48 yield functions are 0.56009, 0.54134, and 0.56488 at the specific plastic work of 3 MPa, and 0.56057, 0.55641, and 0.58394 at the specific plastic work www.nature.com/scientificreports/ of 8 MPa, respectively. It can be seen from the results that the deviation among the yield functions increases with the ratio of shear stress, and the difference is noticeable in the case of σ 12 /σ > 0.5.
In-plane anisotropy of deformation and strength. For polycrystalline metal sheets, the variations of r-value and normalized yield stress ( Y θ ) with uniaxial tensile directions are commonly employed to measure the in-plane anisotropy of deformation and strength. Figure   Nevertheless, it is noticed that the r-values obtained by the virtual tests are a little higher than the experiments; this is might because the present CP model did not consider the potential 111 pencil-glide mechanism of BCC metals 39 , which does not prescribe a specific slip plane for the 111 dislocation glide and thus enables a large number of alternative slip systems. As also stated by Yazawa et al. 40 , limited slip systems increase plastic deformation anisotropy, which is advantageous for realizing higher r-values for BCC steels. In other words, a more advanced physically-based CP might be able to yield an even better prediction of the plastic anisotropy for this type of AHSS. In addition, the experimental data evidence that the variation of r-values with deformation is not noticeable. At the deformation stages with the specific plastic work of 3 MPa and 8 MPa, the r-values are hardly distinguishable in terms of deformation; with the deformation increased to 60 MPa, a small decrease of r can be deduced from Fig. 8, i.e., the decrease of r-values in the directions near DD and the increase of r-values in both RD and TD. The yield functions, as stated, calibrated with the VL generated stress points only, successfully captured the variation tendency of r-values. Especially, the advanced yield function Yld2004-18p predicts the variation of r-values in a rather good agreement with both the virtual tests and experimental results; other two yield functions, i.e., Hill48 and Yld91, slightly underestimate r or the in-plane anisotropy. This is expectable as both have fewer material parameters for describing plastic anisotropy. The Yld91 criterion is a subset version of the  www.nature.com/scientificreports/ Yld2004-18p; the former shares only one-half (either the parameter matrix C or D illustrated in Supplementary Eqs. 16-17 of the parameters of the latter. The quadratic Hill48 yield function also predicts a faithful directionality of the r-values; this is, to some extent, attributed to the plain variation of the r θ curve. Even though the studied DP980 sheet consists of a dual-phase microstructure, it has only two stationary points in the r θ curve; one locates in the RD/TD and another near the DD. The quadratic Hill48 criterion, as anticipated, is capable of capturing this phenomenon. In a word, although only the VL generated stress points were used to identify the anisotropy parameters of the yield functions, all the calibrated functions successfully capture the directionality and variation of r-values. The discrepancies among the predicted results of these functions can also be associated with the intrinsic characteristics of these functions. As shown in Fig. 9, the DP980 sheet exhibits a small strength anisotropy with the maximum variation of the normalized yield stress Y θ below 5%; this can also be demonstrated by the flow stress curves shown in Fig. 5. With the increase of deformation, the strength anisotropy was slightly weakened. Regardless of the rather small variation of Y θ , the stress points obtained by the virtual tests agree well with the experimental data except those in the directions of θ = 30 o and 45 o . All the calibrated yield functions predict a well consistent variation of Y θ with the VL results. It is noted that the Hill48 yield function, which is well-known for its incapacity of capturing the deformation anisotropy and the strength anisotropy simultaneously, gives a satisfactory prediction of the variations of both r-values and normalized yield stresses with tensile direction. Because of the cost and sometimes the restriction  www.nature.com/scientificreports/ of physical experiments, the conventional calibration process of yield functions employs the equal-number experimental data with the parameter number to identify the anisotropic parameters. As stated by Zhang et al. 23 , this routine cannot yield an optimum set of parameters. The Hill48 function, for instance, may predict the deformation anisotropy well but the strength anisotropy poorly when calibrated by the experimental data of r-values, and vice versa. By using enough number (60 here) of stress points with arbitrary loading conditions, the studied yield functions, including the classic quadratic function and advanced non-quadratic functions, calibrated by the full-field CP-based VL, can well predict the deformation and strength anisotropies simultaneously.

Plasticity heterogeneity at the grain level.
To study the role of individual phases in affecting the overall plasticity of the sheet, the microscopic distributions of tensile stress, tensile strain, and r-value were further investigated. Figure 10 presents the contour maps of the deformed RVEs tensile tested in RD at the specific plastic work of 60 MPa; both the as-tested dual-phase RVE (DP-RVE) and the intentionally separated singlephase RVEs (F-RVE for ferrite and M-RVE for martensite) were presented for comparison. The black lines in the DP-RVE depict the phase boundaries. As anticipated, all the RVEs exhibit significantly inhomogeneous distributions of stress, strain, and r-value. The martensite phase, for its higher strength, shows a larger stress level and a smaller strain level in comparison with the ferrite phase. In contrast, the distributions of r-value in the two phases are quite similar and strongly inhomogeneous. The contour maps of the DP-RVE manifest that strain  www.nature.com/scientificreports/ hot-spots occur in some grains' interior and phase boundaries in the ferrite phase, whereas stress hot-spots are mainly in the martensitic grains.
To quantitatively understand the plasticity heterogeneity of the dual-phase microstructure, the histograms of frequency distribution of true strain, true stress, and r-value (shown in Figure 10) were further plotted in Fig. 11. It is noted that we did not distinguish the frequencies of strains above 0.5, stresses above 3000 MPa, and r-values above 5. For example, as shown in Fig. 11a, the large frequency at true strain ( ε ) of 0.25 represents the volume fraction of all material points with ε ≥ 0.25 ; in contrast, the frequencies at ε < 0.25 were calculated using a bin size of 0.005. As shown in Fig. 11a, the strain frequency distribution of the ferrite essentially follows a Gaussian normal distribution with the median of ~ 0.07, whereas that of the martensite is more like a log distribution with the mode of ~ 0.01. It indicates that the ferrite phase, as expected, accommodates most applied deformation and has more homogeneous deformation than the martensite. Li et al. 17 demonstrated that the statistical strain distribution of a single-phase material, either ferrite or martensite, follows the normal distribution. The log distribution of the statistical strain of the martensite in the dual-phase structure is inferred because of the strongly non-uniform deformation partition between the phases. Because of the different strain distributions of the individual phases, the DP-RVE (i.e., the actual DP980 sheet) exhibits a non-probability distribution of strain, suggesting a statically non-uniform deformation in the steel sheet. On the contrary, the stress frequency www.nature.com/scientificreports/ distributions of both phases manifest as the typical Gaussian distribution, as shown in Fig. 11b. This, on the one hand, demonstrates enough number of material points (Fourier grids) and grains (orientations) in the highlyresolved RVE for representing the dual-phase microstructure; on the other hand, it implies the statistically uniform stress distribution of the individual phases. In contrast, the standard deviation of the martensite is obviously larger than that of the ferrite. It implies a more uniform stress distribution in the ferrite. As anticipated, the median stress (~ 1400 MPa) of the martensite is much larger than that (~ 825 MPa) of the ferrite. Besides, unlike the strain frequency distribution, the stress distribution of the DP-RVE is well bell-shaped. The frequency distribution of r-value is rather different from those of strain and stress. For the martensite phase, as shown in Fig. 11c, there are more than 30% material points with r-value close to zero and ~ 13% with r-value above 5.0, and these values for the ferrite is ~ 20% and ~ 11%, respectively. Apart from these two extreme sides, the r-value frequency distributions of the ferrite and martensite phases are greatly similar; both manifest a power-law distribution. These features solidly imply a significant microscopic deformation anisotropy of the DP980 sheet, although its overall normal anisotropy (as shown in Fig. 8) is not very strong. Besides, Fig. 11c displays that the statistical r-value of the martensite is noticeably smaller than that of the ferrite. This phenomenon coincides with the initial crystallographic textures shown in Fig. 4 of the two phases, i.e., the martensite's texture contains a considerable intensity of α-fiber, which weakens r-value of the steel. Besides, the more inhomogeneous deformation in the martensite might also account for its smaller r-value.
It is noted that RVEs tensile tested along other directions exhibit similar results as those presented in Figs. 10 and 11. In summary, these results demonstrate that the CP-based VL with the highly-resolved RVE correlates adequate micro-mechanisms and deformation heterogeneity of the dual-phase microstructure to the macroscopic plastic anisotropy of the DP980 sheet.
Factors affecting the plastic anisotropy of the DP980 sheet. Unlike the FCC structure of the austenitic steel and aluminum alloy, the BCC structure of the ferritic and martensitic steels has no single family of high-density slip planes, which makes BCC metals generally have more complex plastic deformation than the FCC metals. Many researchers reported the good agreement between the predicted and measured r-values for FCC materials such as aluminum and austenitic steel, but poor for BCC carbon steel, ferritic stainless steel, and dual-phase steel [41][42][43] . Compared with the results obtained from mean-field CP simulations 43,44 , the results predicted by the proposed VL show a good agreement with the experimental data, as shown in Figs. 8 and 9. Whereas, it is also noticed that a relatively large discrepancy between the simulated and experimental r-values exists in the DD and TD. To further understand the factors affecting the angular evolution of plastic anisotropy of the dual-phase steel, full-field CP simulations of uniaxial tension along RD, DD, and TD were carried out for single-phase ferritic and martensitic RVEs, respectively; another set of simulations of the dual-phase RVE assuming the same self and latent hardening effect [i.e., q = 1 in Eq. (4)] were also performed. The comparison among the simulation results, in terms of r-values, the average r-value, and the planar anisotropy r , were presented in Figure S2 and Table S2 in the supplementary file.
As shown in Table S2, all the simulations predicted larger r and r than experiments; the original dual-phase RVE yields the most close r and r with experiments, and the RVE with the single martensite phase predicts the largest r-value. It implies that the dual-phase microstructure reduces the plastic anisotropy. As Rollett and Wright 45 and Choi et al. 46 noted, the plastic heterogeneity occurring in a multiphase structure with a soft matrix and hard phases, leads to not only varying magnitudes of strain, but also varying local lattice rotations, all of which tend to disperse texture and reduce anisotropy. It gives one clue to explain the deviation of the simulated r-values shown in Fig. 8: although a full-field CP modeling was adopted in this work, the interaction between www.nature.com/scientificreports/ the two constituent phases, which was enforced by the geometrical constraint among Fourier grids due to deformation compatibility, was much simplified to some extent, especially they have large strength difference and similar volume fraction. On the other hand, the capturing of plastic heterogeneity should be improved by using more advanced CP constitutive theory. A possible way to improve the result is using a physically-based nonlocal crystal plasticity model considering geometrically necessary dislocation enhancement, which can describe a more realistic interaction among heterogeneous microstructures 47 . Besides, the role that crystallography plays at the phase interface has not yet been addressed experimentally and theoretically to much extent. Another factor affecting the simulation results is the latent hardening matrix (or interaction matrix) in Eq. (4). In this work, the latent hardening matrix was referred to the scheme of Peirce et al. 30 , that the hardening ratio q = 1.0 for coplanar slip systems and 1.4 for all non-coplanar slip systems. This scheme is pure phenomenological and simplified 48,49 . As a comparison, simulations with q = 1.0 for all non-coplanar slip systems (the most simplified case) were also carried out. The results presented in Fig. S2 and Table S2 demonstrate that the CP simulations assuming an identical strengthening effect of self-hardening and latent hardening yielded both the larger r and r than the original simulations. With a well-defined interaction matrix considering the different strengthening effects of self-interaction and various collinear and coplanar interactions, one can expect a better performance of CP simulations in capturing the anisotropy of BCC metals. However, the actual latent hardening matrix for a given material is still an open question 50 .
From the aspect of experiments, the mechanical data presented in Figs. 5, 6, 8, and 9 are the macroscopic and bulk mechanical responses of innumerable grains. Although the microstructural RVE was built based upon the EBSD data (as shown in Fig. 4) containing approximately 2000 grains, whereas the EBSD data only represents the crystallographic texture in a local region on the RD-TD plane. Moreover, the thermomechanical operations experienced by the as-received steel generally produce materials with non-uniform spatial distribution of texture. Indeed, texture gradients of both γ and α texture fibers through thickness were commonly reported in cold-rolled and hot-rolled steel sheets 43,51 . As pointed out by Kodukula et al. 43 , one of the key reasons for the large discrepancy between measured and calculated r-values of a BCC steel is the strong texture gradient through thickness. Increasing the ratio of γ − fiber to α− fiber can improve the r-value of BCC steel 52 , and this ratio was reported to decrease from the surface to the center of steel sheets 53 . Figure S3 shows the RD and ND IPFs for the microtextures of studied steel characterized on the RD-TD surface and the thickness center, respectively. It shows that the γ − fiber is significantly intensified in the center. The RVE used in the VL was built without considering texture gradients, and it only incorporated the microtexture data characterized in the center of the steel, which might have the largest intensity of γ − fiber. This could also be a factor resulting in the discrepancy between the simulated and experimental results.

Summary
A virtual laboratory (VL) based on full-field crystal plasticity (CP) modeling was presented to investigate mechanical anisotropy and to predict yield surfaces of multiphase metals. The VL consists of four modules, including a CP constitutive model, a highly-resolved representative volume element of multiphase microstructure, an inverse simulation procedure based on global optimization for identifying the CP parameters of the constituent phases, and a local optimization scheme for calibrating yield functions through a large number of simulated stress points.
Elaborate mechanical experiments and microstructure characterizations were carried out for modeling setup as well as validation of the CP parameters and of the calibrated yield functions. Both the yield loci generated by the VL and predicted by the calibrated yield functions agree well with the experiments, which were enveloped by yield stress points of uniaxial tensions, pure shear, and biaxial tensions with various stress ratios. The deformation and strength anisotropies were finely captured by the VL and the calibrated yield functions. Especially the Hill48 yield function, which is generally incapable of capturing r-value and normalized yield stress ( Y θ ) simultaneously, presents a satisfactory prediction of both r-value and Y θ . This is mainly attributed to the optimized parameters identified from the abundant stress points of arbitrary loading conditions generated by the VL.
Plasticity heterogeneities on grain level were investigated based on the dual-phase RVE, intentionally separated single-phase ferrite and martensite RVEs. It demonstrates that the CP-based VL with highly-resolved RVEs correlates adequately micro-mechanisms and deformation heterogeneity of the multiphase microstructure to the macroscopic plastic anisotropy of the material.
According to the analysis of factors affecting the plastic anisotropy of DP980, the performance of the VL could be improved by introducing a physically-based nonlocal CP model, the role of crystallography plays at the phase interface, a well-defined interaction matrix, and texture gradient through thickness.

Experiments
Microstructural characterizations. The initial microstructure and texture of the as-received sheet were characterized by the SEM electron-backscattered diffraction (EBSD) system, i.e., the VEGA 3 XMU (LaB6) field emission SEM (TESCAN) equipped with an Oxford/Nordlys EBSD detector. The SEM-EBSD system was operated at acceleration voltage of 20 kV and working distance of 20 mm; it scanned the area with size of 100 μm × 100 μm and step size of 0.25 μm, i.e., collected 160,000 data points totally. The data was processed by the AZtec® system (version 2.1, Oxford Inst.) and reproduced with the open-source MATLAB® (version R2018b) toolbox MTEX 38 .

Mechanical tests.
To obtain the parameters of the CP constitutive model and to validate the yield functions determined by the full-field CP-based VL, room-temperature mechanical tests, including uniaxial tension, pure shear, and biaxial tension, were conducted for the as-received sheet. The geometric shape and dimension of these www.nature.com/scientificreports/ specimens are illustrated in Fig. S1a in the Supplementary. The specimens were machined via electrical discharge machining; for the biaxial tensile specimens, seven slits were fabricated on each arm by laser cutting to reduce the geometric constraint on the deformation zone 50 . Uniaxial tensile tests and pure shear tests were carried out quasi-statically on an electronic testing machine (Instron Model 8080 with a load cell of 100kN capacity) equipped with a commercial digital image correlation (DIC) system (ARAMIS), as shown in Fig. S1b, at constant crosshead speeds of 3 mm/min and 2 mm/min, respectively. Uniaxial tensile tests were carried out in the directions of every 15° from RD to TD, and five experiments were carried out for each direction to obtain the repeatable experimental data. The biaxial tensile tests were carried out along four loading paths (with stress ratios σ 11 : σ 22 = 1 : 1, 2 : 1, 1 : 2, and 4 : 3 ) on a biaxial tensile testing machine (MTS BIA5105) as shown in Fig. S1c at a constant equivalent loading rate 0.49 kN/s. Two load cells were employed to measure the tensile force along RD and TD in real-time.
To facilitate the full-field strain measurements of the DIC system, all the specimens were uniformly sprayed with random speckle patterns prior to the tests. The three-dimensional position change of the speckles was photographed by dual high-resolution cameras at a frequency of 2 Hz; the strain field of the deformed specimens then was computed by the commercial DIC software GOM. The photographing was automatically synchronized with the load cell signal by the DIC data acquisition system. www.nature.com/scientificreports/