An informatics method for inferring the hardening exponent of plasticity in polycrystalline metals from surface strain measurements

The investigation of strain hardening in metals is complex, with the outcome depending on experimental conditions, that may involve microstructural history, temperature and loading rate. Hardening is commonly measured, after mechanical processing, through controlled mechanical testing, in ways that either distinguish elastic (stress) from total deformation measurements, or by identifying plastic slip contributions. In this paper, we conjecture that hardening effects can be unraveled through statistical analysis of total strain fluctuations during the evolution sequence of profiles, measured in-situ , through digital image correlation. In particular, we hypothesize that the work hardening exponent is related, through a power-law relationship, to a particular exponent arising from principal component analysis. We demonstrate a scaling analysis for synthetic data produced by widely applicable crystal plasticity models for polycrystalline solids.


Graphical Abstract
When a metal is cold worked, through rolling, drawing or forging, mechanical strength may increase by orders of magnitude, mainly due to the onset of work hardening (Asaro 1983), where dislocations are formed at grain boundaries or interfaces between various coexisting phases.While well exploited, this remarkable physical phenomenon cannot be directly predicted, without case-by-case investigations (Cottrell 2019).For example, the identification of the work hardening exponent requires standardized testing and associated modeling.However, strain hardening critically depends on environmental conditions, such as temperature, pressure, loads, and especially for extreme conditions (Hosford 2010), one needs a large collection of controlled testing that may be impossible (eg. at high temperature or irradiation conditions).In this context, data science and materials informatics (Frydrych et al. 2021) provides tools for the understanding of strain hardening effects with minimal input from experimental set ups.In this paper, we explore the possibility of extracting the strain hardening exponent by the statistical analysis of surface total strain map sequences, that may be extracted through the use of digital image correlation (DIC) methods (Papanikolaou and Alava 2021), that are quite common across scales in the study of metal surfaces (Sutton et al. 2009).We demonstrate, for synthetic data, that features extracted from principal component analysis (Giannetti et al. 2019) of surface strains during uniaxial tensile loading, can be directly associated to the hardening exponent.We propose a functional form that can be tested in experimental settings and in general loading conditions.The characterization of a metal's ability to deform before fracturing, is typically done through various measures of tension tests, such as elongation at fracture, reduction of area at fracture, and strain hardening exponent, typically measured from the engineering yield up to the ultimate yield points.However, all these tests not only depend on the very definition of the yield point (engineering/ultimate), but also necessitate a high level of control and standardization procedure.
DIC methods (Sutton et al. 2009) have emerged as a promising technique for materials informatics applications (Frydrych et al. 2021).DIC uses specimen surface images taken in the initial and deformed configuration to produce displacement and strain fields.DIC is a computationally-intensive technique, but recent progress noticed (Papanikolaou et al. 2019;Mäkinen et al. 2020Mäkinen et al. , 2015;;Yang et al. 2020;Papanikolaou and Alava 2021) that the evolution of DIC total strain profiles has rich information content that may be adequate to infer material properties related to yielding, hardening and eigenstrains, but only a small part of it is utilized for various purposes (Hazeli et al. 2018).Indeed, one may show (Papanikolaou et al. 2019; Papanikolaou and Alava 2021) that the history of deformation of the crystalline sample can be used for a robust definition of the yield point, as the plasticity onset point, in a way that is solely controlled by the data, without engineering conventions.In fact, in Ref. Papanikolaou and Alava (2021), two functional tools were introduced, based on principal component analysis (PCA) and discrete wavelet transforms (DWT), that were shown to capture the onset of crystal plasticity through the analysis of total strain fluctuations during mechanical loading.As a test case, these tools were used on DIC of polycrystalline samples, using synthetically produced data in a phenomenological crystal plasticity model for pure Al.These methods were confirmed to correctly estimate experimental yield stresses in metal alloys (Mäkinen et al. 2022).Beyond yielding, fracture toughness can be also estimated using DIC, typically using artificial neural networks (ANNs) (Cidade et al. 2019;Papanikolaou 2020;Rezaie et al. 2020).Crack detection, measurement, and characterization based on DIC is also done using image processing methods (Gehri et al. 2020) and fatigue crack detection (Strohmann et al. 2021).
In this paper, we focus on polycrystalline metals, analogous to Ref. Papanikolaou and Alava (2021), and small quasi-elastic strains ( < 0.3% ), and the ductile features beyond yielding.We focus on the identification of the hardening exponent n directly out of total strain maps, without the use of local or global stress information.We perform tensile loading (along x) simulations in three-dimensional samples of fixed discretization and polycrystalline structure, with a phenomenological crystal plasticity model, that incorporates well established polycrystalline deformation mechanisms (Roters et al. 2019;Papanikolaou and Alava 2021).We consider synthetic data of uniaxial tensile mechanical deformation in polycrystalline samples with 5µ m local resolution in each directions and grains that are at least 25µ m large in each direction.We vary the microstructural hardening exponent by modifying a particular parameter of the model, the slip interaction parameter h 0 , which ultimately changes the hardening exponent from 0.001 up to 0.3, which is also the typical range seen in experimental data for metals with various degrees of cold work (Silva et al. 2018;Hosford 2010).On all the cases, we implemented the proposed PCA method (Papanikolaou and Alava 2021) by using standard tools in Python (Pedregosa et al. 2011).We show that the plasticity transition is qualitatively modified, but the PCA measures still efficiently track the yield point, and they also allow for the quantification of the hardening exponent n through a constitutive power-law relation.We demonstrate how the constitutive power-law is constructed, and we apply it on the existing data.We conclude by discussing the accuracy of the approach and the possibility of learning additional properties, such as additional hardening information and possibly, damage parameters, but also distinguish various hardening from damage mechanisms.
We study tensile loading in the x-direction, for 3D polycrystalline samples that are periodic in all directions (see Fig. 1(a), and sample dimensions in (x,y,z): (64,64,64) (3D), promoting the perspective of investigating a (0.32mm) 2 sub-mm 3D samples.The crystalline structure of the material is face-centered cubic (FCC) Aluminum (Al), with standard stiffness coefficients (see Table 1, in reference to the cubic coordinates).The examples studied can be readily achieved in modern DIC applications (Sutton et al. 2009) and they follow prior relevant works (Papanikolaou et al. 2019;Roters et al. 2019;Papanikolaou and Alava 2021).The digital surface image collection is assumed to be collected at periodic applied strain intervals (cf.Fig. 1(b)).The samples are periodic in all three dimensions, and for this reason, we utilize all possible y-z surfaces for DIC sampling purposes.Clearly such surface collections cannot be tracked by using typical surface DIC, requiring three dimensional versions, using x-ray tomography.However, this work's impact is solely focused on the demonstration of how surface strain measurements can be connected to realistic crystal plasticity effects.The samples considered in this work display identical polycrystalline texture and similar yield points but distinct hardening exponents n (with the stress σ (ǫ) ∝ (ǫ − ǫ Y ) n ), with the slip-slip interaction coefficient h 0 being altered from 10 6 to 10 12 Pa, while all other sim- ulation parameters remained identical (cf.see also, Fig. 1(c)).Each sample corresponds to a single hardening exponent n and slip-slip interaction coefficient h 0 .
Fig. 1 Synthetic data for phenomenological polycrystal plasticity, hardening effects and strain-map sequences.(a) A 64x64x64 polycrystalline sample is uniaxially loaded under tension along the x-direction, and digital surface image collection is assumed at a 10µm-resolution, assumed to capture a representative volume element, (b) Digital image correlation (DIC) strain maps are assumed to be collected at periodic applied strain intervals on a y-z plane, (c) The samples considered in this work display similar yield points (see Table 1 but distinct hardening exponents n (with the stress σ (ǫ) ∝ (ǫ − ǫ Y ) n ), with the slip-slip interaction coefficient h 0 being altered from 10 6 to 10 12 Pa (see legend for symbols correspondence), while all other simulation parameters remained identical The model (Papanikolaou et al. 2019;Roters et al. 2019) utilizes the phenomenological crystal plasticity theory, capturing slip-based macroscale plasticity, with constitutive laws that are applicable in common metals (Asaro 1983).The model captures in a self-consistent manner the basic physical mechanisms of crystal plasticity, as they take place in most metals, and it captures finite deformations in a cubic grid (Papanikolaou and Alava 2021).All studied samples have 256 grains that are distributed randomly using a Voronoi tesselation in three dimensions (Press et al. 1989) (see also Fig. 2(a) for a y-z surface set of grain orientations) and the model is solved by using a FFT-based method (Papanikolaou et al. 2019).The main evolution of the plastic deformation tensor is governed by: with L p = α γ α s α ⊗ n α , and s, n unit vectors on slip direction and slip plane normal, respectively, while α is the slip system index.Total deformation translates in elastic and plastic ones through F = F e F p .The slip rate γ α is given constitutively through (Asaro 1983;Bronkhorst et al. 1992), (1)  The observed hardening exponent dependence on the slip-slip interaction parameter h 0 for the simulations performed in this work, being apparently non-linear.We acknowledge that the observed non-linearity may be an outcome of the chosen discretization or/and texture, but its study is not related to the subject of the current work where γ0 is the reference shear rate, τ α = S • (s α ⊗ n α ) is the resolved shear stress at a slip resistance g α , with S = [C]E ǫ being the Second Piola-Kirchoff stress tensor, n) in this work (inverse of the strain rate sensitivity exponent m = 1/n ) and g α is the slip resistance for a slip system α .Hardening is provided by the constitutive law (Brown et al. 1989): where h αβ is the hardening matrix: which phenomenologically captures the micromechanical interactions between different slip systems.Here, h 0 = 75MPa, p = 2.25 , and are slip hardening parameters, that are typically assumed to be identical for all FCC systems owing to the underlying dislocation reactions.g β is the resistance to shear on the β slip system, and g β ∞ is the saturated shear resistance on the slip system β (set at g s = 63MPa for all slip systems, see also Table 1), and the shear resistances asymptotically evolve towards saturation.The parameter q αβ is a measure for latent hardening and its value is taken as 1.0 for mutually coplanar slip systems, and 1.4 otherwise, rendering the hardening model anisotropic.
In this work, we utilize the modification of the slip-slip interaction parameter h 0 , as a way to control the hardening exponent in this model.As seen in Fig. 2(b), the hardening exponent changes from 10 −3 to 0.3, as h 0 is varied by 6 orders of magnitude (see also Fig. 1(c) for the changes in the stress-strain curves).We sample the following values of h 0 = 10 6 , 5 × 10 6 , 10 7 , 5 × 10 7 , 10 8 , 5 × 10 8 , 10 9 , 5 × 10 9 , 10 10 , 5 × 10 10 , 10 11 , 5 × 10 11 , 10 12 .Nevertheless, the behavior of n is non-linear with h 0 (see Fig. 2(b)), but it is clearly expected, given the model complexity.For the calculation of n, we utilized the common engineering approach of fitting a linear polynomial function to the logarithm of σ − σ Y , as function of the logarithm of ǫ − ǫ Y , with ǫ Y set to 0.12, and σ Y being set at the value of applied stress at 0.12 loading strain.
Our focus in this work, however, is the use of surface strain maps to predict the change in hardening exponent.For this purpose, we consider the total strain maps of y-z surfaces (perpendicular to the loading direction), at each loading point (seen on Fig. 1(c)), having an image of total strain for every recorded loading step in the simulation, emulating the DIC procedure (Yang et al. 2020;Papanikolaou 2018).The samples are loaded under uniaxial tension along the x-direction.Characteristically, each simulation is composed of ∼ 8000 loading steps, whereas approximately 80 strain maps are considered at 100-step intervals.These conditions resemble the ones that may be generated by DIC for metals under quasi-elastic applied loads (Sutton et al. 2009).
In Ref. Papanikolaou and Alava (2021), it was demonstrated that the identification of the yield point in polycrystals is possible through only total-strain-measurements, (2) p Papanikolaou Journal of Materials Science: Materials Theory (2024) 8:14 without separation of elastic from plastic contributions in total strain profiles, and the model studied was coarse-grained and did not capture various details of active deformation mechanisms (Kubin et al. 1992;Song et al. 2019).The way this feat was accomplished in Ref. Papanikolaou and Alava (2021) was through the application of statistical methods on the raw total strain map sequences across the elastoplastic transition (Papanikolaou et al. 2017).In particular, the use of PCA provided the access to experimentally tractable measures that utilized the Mahalanobis distance, using the basis of N data vectors E k (capturing the norm |F − I| of total distortion), with predominant correlations in the principal directions.
Each E k is a flattened total-strain vector of N = L x * L y , corresponding to an L x × L y total-strain image (see Fig. 3(b)).The success of the method is gauged by the capacity to capture relative fluctuations in just few principal vectors, like the one shown in Fig. 4(b,  c, e, f ).As can be seen in Figs.4(a), PCA works well, in the sense that the cumulative variance recorded in the datasets is captured by the fluctuations seen in at most 4 singular components.
Nevertheless, for near-minimum and near-maximum hardening coefficient ( h 0 = 1 MPa, and h 0 = 10 6 MPa respectively), the plain view of fluctuations seemingly contains little information at total strain 0.3% (see Fig. 3(b, e)).Concurrently, Von Mises stress fluctuations across the sample (seen in Fig. 3(a, d)) display some differences due to the much larger average stress in the sample (see Fig. 1(c)) and the fixed colorbar in the figures.However, the plastic distortion F p is characteristically important for iden- tifying hardening effects, with the plastic distortion maps (|F p − I|) having a clear The Frobenius norm of the plastic distortion (defined as |F p − I| ) surface map for the observed y-z plane at loading strain 0.3% are shown for h 0 = 10 6 Pa in (c) and for h 0 = 10 12 Pa in (f).Notice that yielding is clearly visible in the plastic distortion maps, but not the total distortion ones.Nevertheless, hardening exponent differences, when viewed by the human eye, seem to be isolated on plastic distortion signatures near specific grain boundaries distinction near grain boundaries of yielding grains (compare texture in Fig. 2(a) with plastic distortion maps for h 0 = 1MPa, and h 0 = 10 6 MPa in Fig. 3(c,f ) respectively).
The key conclusion from observing the datasets is that the presence of an always finite, superposing, plasticity-dependent elastic contribution completely masks the effects of the underlying plastic contribution in strain maps, that unveils yielding and hardening effects, keys to plastic localization physics (Asaro 1983).Our findings are also consistent to the understanding that strain-gradients developed near high grain-boundary misorientations lead to emerging stresses and plastic distortions, instead of actual grain orientations (Chen et al. 2010).
Following Ref. Papanikolaou and Alava (2021), these vectors are projected back to the original data vectors through a vector dot product, for datasets that are standardized for the same average value ('0') and standard deviation ('1'): where implies averaging �E k � = 1/N i E (i)  k and i denotes spatial locations, out of N total.
Principal components arise from decomposition of the matrix X = {E k } for k ∈ { recorded loading steps} , with N lines and V columns.The covariance matrix can be diagonalized C = X T X/(n − 1) = V (� 2 /(n − 1))V −1 , and C's V eigenvectors corre- spond to the principal components, while C's 2 eigenvalues are the principal values.
As seen in Fig. 4(b, c) for h 0 = 10 6 Pa (and hardening exponent n ∼ 0 ), and Fig. 4(e, f ) for h 0 = 10 12 Pa (and hardening exponent n ∼ 0.3 ), the first and second principal components are quite different for small and large hardening exponents, showing a strong hardening dependence.However, the differences are not trivial to correlate to hardening exponent values.For this purpose, we consider the post-yield behavior of P ǫ k (2) , which is fitted to a power-law function, following the ansatz: with P m (2) signifying the maximum value of P ǫ (2) , and C 0 , n P2 being fitting constants (see also Fig. 5).
When P ǫ (1) is plotted against P ǫ (2) (cf.Fig. 4(d)), there is an evident parabolic behavior, which is almost universal across all samples (cf.Fig. 6(a)), something that can be noticed when both quantities are plotted with respect to their zero strain value P 0 i .This universal parabolic behavior is further confirmed by plotting the absolute distance of P ǫ (1,2) with respect to their maximum value, showing a quadratic behavior (cf.Fig. 7(a)) |P 2 − P m 2 | ∼ (P 2 − P m 2 ) 2 .In addition, similarly to the results in Ref. Papanikolaou and Alava (2021), there is a standard behavior of P ǫ (1) and P ǫ (2) with respect to the loading strain ǫ , with P ǫ (1) showing a monotonically decaying behavior (cf.(2) shows a behavior similar to order parameter susceptibilities, with a characteristic peak that follows the apparent yield point, as shown in Fig. 6(c).Thus, it appears that the definition of the yield point in (6) Eqs. 7, 8 is uniquely defined and in accordance to what one would estimate through plain view of the stress-strain behavior.
Ultimately, one can consider the behavior of n P2 with respect to the hardening expo- nent n, which can be seen in Fig. 7(b) for all samples studied, which shows a power law behavior.We conjecture that there is a general constitutive connection between the exponents n P2 and n: with δ ≃ 0.5 (see line on Fig. 7(b)) for the material class studied in this work and C 1 , C 2 being fitting constants.
In summary, we demonstrated that the identification of the hardening exponent in polycrystals is possible through only total-strain surface-measurements, without the need to separate elastic from plastic contributions from total strain profiles.The current model study serves only to provide a demonstration and is not destined to be directly comparable to experimental setups, a topic of interest in a future work (Mammadli et al. 2024).In the performed simulations, each pixel's linear dimension can be thought of corresponding to 5µ m (the scale of a representative volume element in the microstructure) and the linear size of tracked areas approaches 0.35mm.We find that predictability of ( 10) In fact, The model studied is coarse-grained and does not capture detailed deformation mechanisms (Kubin et al. 1992;Song et al. 2019).An interesting future direction involves the extension of this method to creep and fatigue conditions, that can thus have industrial relevance.In addition, one may find analogous dependencies for damaged specimen synthetic data.It is important to identify the dependencies, and relevant scaling functions (Papanikolaou et al. 2017), both being the subject of future works.

Fig. 2
Fig. 2 Polycrystallinity and hardening exponent variability.a Grain orientations on the observed y-z surface, perpendicular to the loading-axis x (given the representation through Bunge Euler angles with respect to the x-axis projection is cosφ 1 cosφ 2 -sinφ 1 sinφ 2 cos ), and the RGB color index being proportional to (φ 1 , φ 2 , �) .bThe observed hardening exponent dependence on the slip-slip interaction parameter h 0 for the simulations performed in this work, being apparently non-linear.We acknowledge that the observed non-linearity may be an outcome of the chosen discretization or/and texture, but its study is not related to the subject of the current work

Fig. 3
Fig. 3 Polycrystalline texture and mechanical response.Samples studied in this work are on 64x64x64 (3D) grids, with 256 grains in each microstructure.The Von Mises stress surface map for the observed y-z plane (having the microstructure shown in Fig. 2(a)) at loading strain 0.3% are shown for h 0 = 10 6 Pa in (a) and for h 0 = 10 12 Pa in (d).The Frobenius norm of the total distortion (defined as |F − I| ) surface map for the observed y-z plane at loading strain 0.3% are shown for h 0 = 10 6 Pa in (b) and for h 0 = 10 12 Pa in (e).The Frobenius norm of the plastic distortion (defined as |F p − I| ) surface map for the observed y-z plane at loading strain 0.3% are shown for h 0 = 10 6 Pa in (c) and for h 0 = 10 12 Pa in (f).Notice that yielding is clearly visible in the plastic distortion maps, but not the total distortion ones.Nevertheless, hardening exponent differences, when viewed by the human eye, seem to be isolated on plastic distortion signatures near specific grain boundaries

Fig. 5
Fig. 5 Method for estimating the hardening exponent by using DIC data: An example.The key aspect of the proposed method in this work is the post-yield fit of the decay of the P 2 observable, solely calculated from experimental surface strain maps.Through the fit, one may estimate the exponent n P2 , which captures the decay ( P ǫ 2 = P m 2 − c ǫ − ǫ PCA Y

Table 1
Model parameters chosen in this work.The key parameter h 0 is varied for altering the hardening exponent