Combining Single-and PolyCrystalline Measurements for Identification of Crystal Plasticity Parameters : Application to Austenitic Stainless Steel

Crystal plasticity finite element models have been extensively used to simulate various aspects of polycrystalline deformations. A common weakness of practically all models lies in a relatively large number of constitutive modeling parameters that, in principle, would require dedicated measurements on proper length scales in order to perform reliable model calibration. It is important to realize that the obtained data at different scales should be properly accounted for in the models. In this work, a two-scale calibration procedure is proposed to identify (conventional) crystal plasticity model parameters on a grain scale from tensile test experiments performed on both single crystals and polycrystals. The need for proper adjustment of the polycrystalline tensile data is emphasized and demonstrated by subtracting the length scale effect, originating due to grain boundary strengthening, following the Hall–Petch relation. A small but representative volume element model of the microstructure is identified for fast and reliable identification of modeling parameters. Finally, a simple hardening model upgrade is proposed to incorporate the grain size effects in conventional crystal plasticity. The calibration strategy is demonstrated on tensile test measurements on 316L austenitic stainless steel obtained from the literature.


Introduction
Crystal plasticity theory has been extensively employed in the research of metal plasticity.Since the introduction of the concept of dislocations [1,2], many advances in the understanding of slip systems, work hardening and texture evolution have been made.These findings led to the development of crystal plasticity (CP) models and their implementation into commercial finite element (FE) solvers.Recent developments in crystal plasticity finite element modeling (CPFEM) have enabled the investigation of a variety of aspects concerning the heterogeneous plastic response [3][4][5][6][7] and texture evolution [5,[8][9][10][11][12][13][14] of polycrystalline metals under general boundary conditions.
To accurately predict local stresses and strains with CPFEM simulations, proper identification and reliable calibration of the underlying CP models are obviously needed.For example, there is a growing tendency of building [15] and using FE models that resemble true grain topology, true texture and/or realistic loading conditions with the intention to compare the predictions of emerging local phenomena of the underlying CP models (e.g., grain boundary stress concentrations [16][17][18], plastic strain localization and the formation of slip bands [19,20], crack initiation and propagation [6,19,21]) directly with measurements.In such cases, appropriately-calibrated constitutive laws would certainly be required, which, however, is not a trivial task.Due to the potentially large number of constitutive modeling parameters (typically ∼10 in conventional CP models) and CPU time-consuming CPFEM simulations (especially true for large polycrystalline models), the identification strategy should be cleverly designed and automatized.
In principle, CP parameters should be identified upon measurements performed on a (length) scale of the underlying constitutive model: if CP behavior is modeled on a grain scale, experiments on single crystalline specimens should be considered in the calibration.In this way, the effects of deformation mechanisms not considered in the CP model would be automatically avoided from the data.However, the availability of measurements on a required (length) scale may be limited, and usage of experimental data on other (usually larger length) scales may become tempting or even necessary.It is important to realize that such data should be properly adjusted before they can be used in CP model calibration.
To overcome the issue of using data on different length scales, novel approaches were recently developed to determine reliable single crystal material parameters by using high-energy synchrotron X-ray diffraction measurements to measure grain scale deformation behavior within a polycrystalline aggregate during in situ mechanical loading [22][23][24].The diffraction data obtained directly on a grain scale were used together with crystal-based finite element simulations to identify single crystal elastic moduli [22,23] and slip system strength evolution [24] in titanium alloys.
In the present work, a reliable calibration strategy is proposed for the identification of CP model parameters when experimental tensile data are combined from two length scales: single crystals (meso scale) and polycrystals (macro scale).It is argued that grain boundary strengthening effects due to finite grain sizes should be effectively removed from the polycrystalline data in order to be considered in CP model calibration simultaneously with single crystalline data.In this respect, a well-known Hall-Petch relation [25,26], is used to correlate flow stress σ with average grain size D .Both σ 0 ( ) and K( ) are strain-(and temperature-) dependent macroscopic material properties; at yield stress (σ = σ y ), σ 0 is usually interpreted as the frictional stress opposing the glide of dislocations and K a measure of grain boundary resistance to the transmission of slip from one grain to the next.The proposed calibration strategy is applied to the conventional, length-scale independent CP model, originally proposed by Hill and Rice [27], with a more advanced constitutive hardening law, initially proposed by Bassani and Wu [28] for Stage I and II hardening and later upgraded by Bassani [29] to include also Stage III hardening.An automated fitting approach is developed where model parameters are limited to discrete values with prescribed precision in order to speed up the calibration procedure.The feasibility of the approach is furthermore expanded by identifying an FE model of the representative volume element (RVE) of the microstructure, which is computationally inexpensive and, in this respect, comparable to single crystalline FE models.The identification of the parameters is performed with respect to tensile test measurements on AISI 316L austenitic stainless steel obtained from the literature: single-crystal measurements at 22 • C for the three tensile directions [001], [111] and [123] were published in [30] and polycrystal measurements for different average grain sizes at 24 • C in [31].Finally, a simple upgrade of the hardening model is proposed to include the length scale in the conventional CP model.

Constitutive Model
In the present work, conventional single crystal plasticity [27] is employed with anisotropic elasticity governed by Hooke's law.The corresponding elastic constants for 316L stainless steel are considered [13,32]: c 11 = 204,600 MPa, c 12 = 137,700 MPa, c 14 = 126,200 MPa.The model assumes that plastic flow in single crystals takes place through slip on prescribed slip systems.The FCC crystal considered here has 12 slip systems available for plastic deformation where each slip system α is defined by its slip plane normal n α i and slip direction m α i .It is assumed that deformation takes place through dislocation glide, and the evolution of the plastic velocity gradient L p ij is given by: where the summation is performed over all active slip systems α and γα represents shear strain rate.The crystalline slip is assumed to obey Schmid's law where γα is assumed to depend on the current Cauchy stress σ ij solely through the Schmid stress τ α defined as: The shear strain rate for a visco-plastic flow can be expressed in the form of a power law [33]: which is a reasonable approximation for metals at room temperature, ordinary strain rates and pressures [34].Here, γα 0 represents a reference strain rate; n is the rate sensitivity index; and s α is the current strength of slip system α.Both γα 0 and n were set to constant values of, respectively, 10 −4 s −1 and 20 to reproduce negligible viscous behavior.The strength rate can be expressed in the following general form: where h αβ is a slip hardening matrix describing the hardening of the slip system α due to slip activity on any other slip system β.
In this work, an empirical constitutive hardening law of Bassani and Wu [28,29] is employed with the following form for hardening moduli (see Figure 1a for the graphical interpretation), Here, h 0 and h s define the hardening slopes immediately following initial yield and during easy glide within Stage I hardening, respectively.Parameter τ 0 is yield stress, equal to the initial value of the strength, s α (0) = τ 0 , and τ s is the saturation stress or a reference stress where large plastic flow initiates.The off-diagonal terms h αβ for α = β are defined through the small parameter q.The second term in the expression for h αα deals implicitly with cross-hardening that occurs between slip systems during Stage II.Here, f αβ is an interaction matrix that depends on the nature of the junctions formed between slip systems α and β, and γ 0 represents the amount of slip after which the interaction between slip systems reaches the peak strength.In general, the interaction matrix f αβ for FCC crystals involves five different parameters whose values depend on the relative strengths of the junctions [28].However, in this study, the matrix has been simplified to a single value, f αβ = f 0 .
The hardening law in Equations ( 6) and ( 7) includes the effects of Stage I and Stage II hardening.Bassani [29] later proposed a way to include the effects of Stage III hardening by assuming that h s depends on total accumulated slip Γ = ∑ α | γα |dt on all systems, where γ I I I 0 is approximately the accumulated slip at the onset of Stage III deformation, while h I s and h I I I s define the hardening slopes at the beginning of Stage I and Stage III deformation, respectively.The extended version of the law, Equations ( 6)- (8), was used throughout this work.The hardening constitutive equations were implemented numerically into finite element solver Abaqus [35] by upgrading Huang's user material subroutine (UMAT) [36] on the parts of the code that address self hardening.The CP model was implemented for finite deformations and rotations using the forward gradient time integration scheme and linear incremental formulation.The general static implicit method was used in Abaqus to solve the constitutive equations.

Two-Scale Calibration Strategy
The proposed strategy for the two-scale calibration of the CP model comprises the following main assumptions and ideas:

•
The hardening parameters are identified on a grain scale by fitting the calculated tensile responses to corresponding experimental stress-strain curves.Simultaneous fitting is performed against single and polycrystalline tensile data.

•
Before being used in the calibration, the length scale effects due to grain boundary strengthening need to be subtracted from the measured polycrystalline data.The Hall-Petch relation is used to obtain the adjusted (infinite-grain-size) macroscopic tensile stress (σ 0 ( ) in Equation ( 1)).

•
The RVE model of the polycrystalline microstructure is identified and used in the identification procedure.To allow for fast simulations during calibration iterations, a small representative polycrystalline model is determined, containing a comparable number of finite elements as single crystalline FE models.

•
The automated optimization procedure for hardening parameters is introduced where parameters are constrained to take only discrete values with a prescribed precision (step size).

Subtracting Grain Boundary Strengthening Effects from Raw Polycrystalline Data
The adjusted macroscopic tensile stress curve, free of grain boundary strengthening effects, was extracted from the published data for polycrystalline 316L stainless steel [31].In [31], tensile specimens of 316L stainless steel with average grain sizes D in the range of 3.1-86.7 µm were deformed in tension to 0.34 strain at temperatures of 24, 400 and 700 • C and a strain rate of 10 −4 s −1 to investigate the Hall-Petch relationship; Equation (1).The experimental stress-strain curves obtained at 24 • C are shown in Figure 1b, where in general, higher flow stresses are observed for smaller grain sizes.After performing Hall-Petch type plots between σ( ) and D −0.5 at several constant strain levels at 24 • C, it was found that up to ∼0.05 strain, the Hall-Petch relationship exhibits a bi-linearity with larger slope K and smaller friction stress σ 0 for fine grains, D ≤ 7.3 µm, than for coarse grains (please refer to Figure 4a in [31]).For strains larger than ∼0.05, a single Hall-Petch slope was observed.
For the purposes of this study, stress σ 0 ( ) was extrapolated from tensile data at 24 • C shown in Figure 1b.The extrapolation was performed at seven macroscopic strain levels: 0.002, 0.01, 0.02, 0.05, 0.10, 0.24 and 0.34.For ≤ 0.05, σ 0 ( ) was obtained by fitting a Hall-Petch relation over the coarse grain size range, D > 7.3 µm, while the entire grain size range was used to obtain σ 0 ( ) for ≥ 0.10.The extrapolated stress σ 0 ( ) for 316L stainless steel is shown by the red line in Figure 1b, representing the limit D → ∞.

Representative Volume Element
To identify an RVE for a 316L stainless steel specimen of gauge size 12.5 × 5.4 × 25.4 mm 3 [31], several polycrystalline FE models with regular partitioning of the modeling space were built.Aggregates with different numbers of grains N g , different mesh densities, denoted by the number of elements per grain N eg , as well as different random initial grain orientations were analyzed in order to identify later the corresponding continuum limit, N g → ∞ and N eg → ∞, of the macroscopic tensile response.Table 1 lists all of the aggregate models used in the study, and Figure 2 shows a few examples.To allow for systematic study, the variability of grain shapes was intentionally omitted by assigning a cubic shape to all of the grains.It seems reasonable to assume that such a rather unrealistic grain topology should not affect substantially the calculated macroscopic tensile response.In fact, this is justified later by performing comparisons with more realistic Voronoi tessellations.The grains were meshed by linear hexahedral elements (named C3D8 in Abaqus [35]) using finite element solver Abaqus.The cubic shape of the elements (as well as grains) was preserved at all mesh refinings.Regarding boundary conditions, an incremental tensile displacement was applied along the axial direction (Z axis) to all of the nodes on the front surface, while the nodes on the back surface were constrained to have zero axial displacement.The applied strain was set to 32% and strain rate to 10 −4 s −1 , in accordance with [31].Using general static implicit scheme in Abaqus, the applied strain and strain rate were implemented by setting total tensile displacement and assigning total simulation time, respectively.All four side surfaces of the model were assumed free of the constraints.
In Abaqus, the grains were defined through sections of common crystallographic orientations.Random initial grain orientations were assigned to all of the models.Models with the same number of grains shared the same set of grain orientations.In aggregates with a smaller number of grains (N g = 64 and 512), however, twenty different random set orientations were generated to estimate an average tensile response.A random set that provided the closest response to the average one was then selected to be used in the extrapolation analysis.As expected, models with more grains (N g ≥ 4096) showed practically no variability with respect to different initial random sets.
The continuum limit of the realistic 316L stainless steel specimen was identified through the extrapolation of calculated tensile curves on various polycrystalline aggregates with different numbers of grains and elements per grain.Figure 3 presents the simulated stress-strain curves of some of the models from Table 1 calculated with the same set of CP parameters.A large scatter of calculated tensile curves is observed, however, with the following general tendencies: (i) a stiffer tensile response is obtained for increasing number of grains N g (at constant N eg ), and (ii) a softer tensile response is obtained for increasing number of elements per grain N eg (at constant N g ).The latter behavior (ii) can be understood in terms of putting more degrees of freedom toward the grains, thus effectively making them softer.As linear elements, in general, assume a constant strain (and stress) within the element, possibly large local strain gradients, physically present due to complex grain loading from the surrounding (anisotropic) neighborhood, may not be accurately resolved in the simulation when using too small a number of elements per characteristic length (i.e., N eg ).As a consequence, the whole grain may appear as too stiff (unable to deform adequately to prescribed loads, thus providing, on average, higher stress).While Behavior (ii) is a numerical feature, Behavior (i) is related to the physics of the problem and can be connected to the probability of finding a cross-section along the specimen with the softest grains; with smaller N g , soft cross-sections are more likely to be found, thus providing lower tensile response.Although similar in effect, this behavior should not be misinterpreted as a grain boundary strengthening effect, which was avoided by design by using length-size independent CP law.
Figure 3. Calculated tensile responses of some of the aggregate models from Table 1 with (a) fixed N g = 64 and increasing N eg and (b) fixed N eq = 1 and increasing N g .The two arrows depict the values of strains for which a convergence analysis was performed.
In Figure 4, the convergence analysis of the calculated tensile stresses from Figure 3 is performed at true strains 33 = 0.1 and 0.2.The goal of the analysis is to identify a converged stress value by extrapolating the calculated stresses to the limiting case of N g → ∞ and N eg → ∞.In Figure 4a,b, the two opposite stress tendencies for increasing N g and N eg are clearly visible.However, they both seem to nicely follow a power law behavior.For this reason, a power law fitting function was constructed to account simultaneously for both variables N g and N eg , where σ * 33 and a i are fitting parameters to be adjusted to the stress data.The resulting fits (see fitted dashed lines) show an excellent agreement between the calculated data points and Equation (9).The results of the fits are summarized in Table 2.
Table 2. Results of the fitting of Equation ( 9) to the calculated stress data of Figure 3  The values of σ * 33 were finally used to identify the RVE model from the list in Table 1 that provides the closest stress value to σ * 33 for all considered strains.According to the extrapolations performed at strains 0.1 and 0.2, the RVE model is selected to be a model with N g = 64 grains and N eg = 8 elements per grain.From Figure 4, it can be estimated that stresses of the proposed model should not deviate more than 10 MPa from the stress σ * 33 .This model is also small enough (N e = 512 of total elements) to allow short simulation times in the calibration procedure.

Single-Crystal Models
Single-crystal models were built from 48 (3 × 2 × 8) quadratic hexahedral elements with reduced integration (named C3D20R in Abaqus [35]) to model the realistic specimen gauge section of 3 × 1.5 × 8 mm 3 of 316L stainless steel used in [30].Three directions of straining, [001], [111] and [123], were modeled by assigning appropriate material orientations.The finite element mesh was checked for convergence for all three tensile directions.In this way, an estimated accuracy of a few MPa in macroscopic tensile response was achieved with a 48-element model for the [123] direction of straining.Tension tests were simulated with the displacement control analysis using a fixed strain rate of ˙ = 5 × 10 −5 s −1 according to [30].The displacement boundary condition was specified in the axial direction (Z axis) for all of the nodes on the front surface, and the nodes on the back surface were constrained to have zero axial displacement.For the [123] direction of straining, the global rotation of the model around tensile axis was prevented by additionally constraining two pairs of edges on the front and back surfaces.

Automated Optimization of Model Parameters Defined on Discrete Range
The automated optimization procedure was originally proposed in [37] and later used in [13] for the calibration of the simpler Stage I and II hardening law of Bassani and Wu [28] against single crystal measurements on 316L stainless steel.In the present work, the procedure is upgraded to account also for polycrystalline data and optimized for speed by using a discrete parameter set.
In the optimization procedure, the CP parameters are identified by fitting simultaneously the three calculated stress-strain curves (strained along [001], [111] and [123] directions) and the fourth, polycrystalline tensile curve, to the corresponding single-and adjusted poly-crystal measurements on 316L stainless steel.The optimal values for the parameters are obtained by minimizing the χ 2 merit function, where σ ori,i and ori,i are, respectively, the measured true stress and true inelastic strain for a given direction of straining denoted by ori.The σ ori ( ori,i , P) is the calculated true stress component along tensile direction (Z axis) obtained for a given crystal orientation (or polycrystal) and parameter set P.
In addition, a flexible fitting domain is introduced by setting a maximum strain in the calculation of χ 2 (addressed by variable m).Here, four fitting domains were used with corresponding maximum nominal strains e max = 0.1, 0.2, 0.3 and 0.4.Minimization of χ 2 with respect to model parameters p i ∈ P was performed using Powell's minimization method [38].In this respect, a two-way communication between Abaqus and the minimization subroutine was established through a separate program, which enabled an automatic identification of the parameters.To speed up the iteration procedure, the parameters p i were furthermore constrained to take only discrete values p i + j∆p i on a given range p i,min ≤ p i ≤ p i,max .Here, a fixed step was selected to be 2% of the initial parameter value, ∆p i = 0.02p (0) i .Nine fitting parameters of the hardening law are summarized in Table 3.
Typically, 10 3 FE simulations were needed in Abaqus to achieve the convergence of Powell's method.For a given mesh density (i.e., 48 elements for the single-and 512 elements for the poly-crystal model), this resulted in 1-2 days of CPU time.Several runs with different initial parameter sets were tested to approach the global minimum of χ 2 .

Results of Model Calibration
Results of the two-scale calibration procedure are summarized in Figure 5 for four fitting domains studied in this work.In general, very good agreement is observed between the fitted lines and corresponding measurements in all considered domains.This strongly suggests that the proposed CP model with the Bassani and Wu hardening law [29], denoted as BW3, proves to be adequate to describe the tensile behavior of 316L stainless steel at room temperature both at single-and poly-crystalline scales.
Naturally, not all domains provide equally good fits.For example, in the shortest domain with e max = 0.1 (Figure 5a), an almost perfect fit is obtained with the average stress deviation of only χ 2 = 6.7 MPa and the corresponding parameter set P 1 given in Table 4.However, extending the curves beyond e max = 0.1 (dashed lines in Figure 5a), but computed with the same set P 1 , produces poorer results especially in [111] direction.This rather large misfit in the [111] direction at > 0.1 is successfully suppressed when using larger fitting domains, e max ≥ 0.2.However, few other discrepancies emerge on the other two tensile directions ([001] and [123]).It is interesting to note that the agreement with the (adjusted) polycrystalline data remains very good up to e max = 0.4.Since no information on the measurement error is given neither in [30], nor [31], the obtained fits for e max = 0.2, 0.3 and 0.4 are still considered as very good (note that χ 2 < 21 MPa), thus successfully validating the employed BW3 model.The corresponding parameter sets P 2 , P 3 , P 4 can be found in Table 4.
Small variations among the parameters of the three sets (P 2 , P 3 , P 4 ) suggest that the parameters are well determined and thus belong to the same "material state".The largest uncertainty can be assigned to parameter h I s , which shows variations from 71 MPa in P 2 to 125 MPa in P 4 .This may indicate that, in addition, other tensile directions would be required in the fitting procedure in order to reduce the ambiguity of h I s .Moreover, a relatively small value of γ I 0 ∼ 0.02 and a high value of f 0 ∼ 0.67 indicate a short Stage I region and hard activation of the cross slip during Stage II hardening.This agrees with the behavior of low stacking-fault energy steels [39].A low value of q ∼ 0.2 also suggests that self hardening dominates over latent hardening in 316L stainless steel at room temperature.4.

Introduction of Length Scale into Crystal Plasticity
As conventional crystal plasticity is length-scale independent, here, a simple modification of the hardening law is suggested to incorporate grain boundary effects in a polycrystalline simulations following the Hall-Petch relation; Equation (1).The idea is to introduce an empirical rather physics-based modification that will produce a desired macroscopic polycrystalline response [41][42][43][44].
It is assumed that hardening moduli h αβ remain length-scale independent, and the inclusion of Hall-Petch relation is carried out by adding a new term to Equation ( 5) for the strength rate and initial strength of the material, The introduced grain boundary strengthening is assumed to be isotropic (i.e., slip system independent): with k(Γ, D ) being a grain-scale Hall-Petch parameter that is assumed to depend on total accumulated slip on all systems, Γ = ∑ α | γα |dt, and also on average grain size D .This last dependence on D is introduced to allow for a non-standard Hall-Petch relationship (e.g., bi-linearity observed in 316L stainless steel at room temperature [31], see also Section 3.1).For the standard Hall-Petch relation, k = k(Γ).
In the limit where k is weakly dependent on Γ, ṡHP ∼ 0, and the proposed modification merely shifts the macroscopic tensile curve σ 0 ( ) in the vertical direction, which is in accordance with experimental observations on 316L stainless steel at least for 0.1 and D 7.3 µm (see Figure 1b).This suggests that k may be accurately expressed in terms of the powers of Γ with only a few lowest terms retained, e.g., k = ∑ j k j Γ j ≈ k 0 + k 1 Γ.In this way, the modified hardening law, Equations (11) and (12), can be simplified to: Note that this linearized approximation in Γ is similar in effect to the law introduced in [43,44].Parameters k 0 and k 1 , which may also depend on D , can be identified by (second) fitting of the simulated polycrystal response to a given experimental curve with finite D ; see Figure 6a.The remaining parameters of the BW3 model used in Figure 6a were taken from set P 4 of Table 4. Very good agreement is observed with all tensile measurements on 316L stainless steel [31] on the whole strain domain, which confirms the validity of the proposed modification in Equations ( 14) and (15).However, the obtained k 0 -and k 1 -values, optimized separately for each D , are not unique, but scattered with corresponding averages and standard uncertainties: k 0 = (79 ± 16) MPa √ µm and Part of this scatter can be attributed to the bi-linear Hall-Petch relationship observed in this material for 0.1.It is expected that the scatter can be further reduced if tensile measurement errors (not reported in [31]) would be also accounted in the fits.
To avoid the additional fitting performed in Figure 6a, k(Γ, D ) may be in general related directly to experimental K( , D ) from Equation (1) as: by using grain-averaged Taylor factor M for a random (untextured) polycrystal under tension: Here, σ 0y represents the macroscopic yield stress of a polycrystal with D → ∞, and p is a macroscopic plastic strain (note, however, that p ≈ ).In the limit of homogeneous strain (Taylor model), M ∼ 3 for FCC polycrystals, while M ∼ 200 MPa/82.8MPa ∼ 2.4 can be estimated from the calculated data for 316L stainless steel.The value M ∼ 2.4 should in fact apply to any random FCC polycrystal with slip-driven plastic flow.
The modified hardening law, Equations ( 11) and ( 12), can finally be expressed solely in terms of experimentally determined Hall-Petch parameter K ≡ K( , D ) as: For demonstration purposes, K(0, D ) and ∂K/∂ were estimated from [31] for 316L stainless steel at room temperature only for larger grain sizes, D ≥ 7.3 µm, to avoid the bi-linearity of the Hall-Petch relation: K(0, D ) ∼ 135 MPa √ µm, ∂K/∂ ∼ 0 for ≤ 0.05 and ∂K/∂ ∼ 750 MPa √ µm for > 0.05.These values were used in polycrystalline simulations along with the parameters of set P 4 of Table 4 to calculate the tensile responses shown in Figure 6b.Although no additional fitting was required in the calculations, still a very good agreement is observed with corresponding tensile measurements, thus confirming the validity of the proposed modification in Equations ( 18) and (19).It is to be noted, however, that the introduced length scale in the hardening law is intended to be used when performing polycrystalline simulations with CP parameters determined on a single crystal scale (i.e., without grain boundary strengthening effects).In the calibration itself, as proposed in Section 3, the length scale effects are excluded from the experimental data; therefore, no modification of the hardening law is needed.14) and ( 15), calculated with newly identified parameters k 0 and k 1 (shown in legend in units of MPa √ µm) and previously identified set P 4 from Table 4.A comparison is shown with tensile measurements on 316L stainless steel at room temperature [31].(b) Same as (a), but using Equations ( 18) and ( 19) for the modified BW3 hardening model and experimental Hall-Petch slope K( , D ) (shown in the legend in units of MPa √ µm) from [31].

Comparison with Simpler Hardening Models
The presented results can be put in the perspective of previous studies on 316L stainless steel.In [37], the original hardening model of Bassani and Wu with Stages I and II [28] (BW2) was calibrated with respect to single crystalline measurements alone.To allow for a comparison, a similar study was conducted here using the same BW2 model, Equations ( 6) and (7) with h s = h I s and γ 0 = γ I 0 , but with the proposed two-scale calibration approach to account also for polycrystalline data [31].Table 4 shows four best-fitting parameter sets (P 5 , P 6 , P 7 , P 8 ) obtained on four strain ranges e max .The sets agree to a high extent (within ∼10%) to corresponding sets reported in [37], where, however, a larger variability of the parameters was observed across different fitting domains.This can be readily explained by noting that fitting in [37] was performed with three (and not four) tensile curves.The fact that the inclusion of the additional polycrystalline fit reduces the ambiguity of the parameters without producing additional off-set is in favor of the proposed two-scale calibration strategy.
Furthermore, a comparison between the BW2 and BW3 fitting sets in Table 4 reveals high similarity also among the two hardening models: P i ∼ P i+4 for all i = 1, 2, 3, 4.Although with a slightly bigger (∼1 MPa) average stress deviation χ 2 , the BW2 tensile curves (shown in Figure 7a for e max = 0.2) are practically the same as those presented in Figure 5 for BW3.It is interesting to note that in BW3 fits, the optimal h I I I s is similar in amplitude to h I s , which, according to Equation ( 8), translates BW3 close to the BW2 limit.This may indicate a (very) slow Stage II to Stage III transition with slip in 316L stainless steel at room temperature.Therefore, at these conditions, both BW2 and BW3 hardening models are shown to provide an equally good description of material deformation.
Table 4 contains also the calibrating sets (P 9 , P 10 , P 11 and P 12 ) of the relatively simple hardening model of Peirce, Asaro and Needleman [40] (PAN) with the following hardening moduli, The parameters of PAN have an analogous meaning as in BW2 and BW3, and Γ = ∑ α | γα |dt denotes the total accumulated slip on all systems.Although a reasonably good fit is obtained for e max = 0.1 ( χ 2 = 9.9 MPa), the results for e max ≥ 0.2 are already less convincing ( χ 2 ≥ 27.5 MPa).The corresponding tensile lines for e max = 0.2 are shown in Figure 7b 4 and 5.
The importance of removing grain boundary strengthening effects from polycrystalline data before using it in the two-scale calibration procedure is quantified in Figure 7c,d and Table 5.In this respect, the adjusted polycrystalline tensile curve σ 0 ( ) was replaced by raw measurements obtained on the 316L stainless steel specimen with the average grain size D = 3.1 and 86.7 µm [31].The raw tensile data were then used simultaneously with the three single-crystalline curves in the calibration of the BW3 model for e max = 0.2.The resulting best-fit parameter set P 14 for D = 86.7 µm is very similar to P 2 obtained on the extrapolated D → ∞ model (Table 5).Furthermore, the corresponding tensile curves in Figure 7d are shown to agree well with the measurements.This indicates that grain boundary strengthening effects in 316L stainless steel with D = 86.7 µm are relatively weak and may as well be ignored at room temperature.However, the opposite is observed in Figure 7c for D = 3.1 µm.Here, the grain boundary strengthening effects are too strong to be ignored in the hardening model; therefore, a poor fit is realized.The corresponding parameter set P 13 is obviously not representative of 316L stainless steel.As demonstrated above, the polycrystalline tensile data need to be properly adjusted for reliable two-scale calibration of the length-scale independent hardening model defined on the grain scale.However, it has to be noted that raw polycrystalline data may, of course, be used directly in a single-(macro) scale calibration.In that case, high uncertainty of the fitting parameters is expected given that the number of modeling parameters is large enough (e.g., ∼10 as in BW2) and the fact that they are obtained by fitting to one single macroscopic curve (in the case of untextured material).Table 6 lists a few such examples for 316L stainless steel taken from the literature [6,14,15,32,45].Indeed, a large scatter of the hardening parameters is observed when compared to the results of Table 4.It is to be noted that model parameters from Table 6 fail to reproduce the single-crystal tensile responses of 316L stainless steel.Moreover, it is easy to anticipate that large uncertainty of modeling parameters may also induce large uncertainty of predicted local fields calculated within a polycrystalline model (see the next section for more details).Anyhow, finding a "true" parameter set from raw polycrystalline data is certainly possible, but may require additional validation with other complementary experiments [46].

Local Fields
In the following, the uncertainty of local stress fields is estimated in a more realistic Voronoi polycrystal FE model with large grain sizes (i.e., diminishing grain boundary strengthening effects) when using different hardening models with corresponding best-fit parameters P 2 , P 6 and P 10 from Table 4.The goal is to quantify the differences in the local von Mises stress of a polycrystal when using different hardening models calibrated to the same material, as shown in Figures 5b and 7a,b.
In this respect, the RVE model from Figure 2a with 64 grains and 512 elements was replaced with a Voronoi aggregate model with 1000 grains and 1.4 million quadratic tetrahedral elements (named C3D10 in Abaqus [35]).Software package Neper [47] was used to build the model geometry and the corresponding FE mesh.The FE model is shown in Figure 8a.The average grain in the model was assumed isotropic in shape (i.e., not elongated in any direction).All of the grains were assigned the same material properties, but random crystallographic orientations, thus providing zero texture to the aggregate.Grain boundaries were not modeled explicitly.In fact, various other Voronoi models were constructed with different numbers of grains and mesh densities to assure that the tensile response of the presented model with 1000 grains is a good approximation of the continuum limit.In Figure 8b, the tensile responses of the Voronoi model are shown for the three calibrated hardening models with corresponding sets P 2 , P 6 and P 10 from Table 4.As already mentioned in a previous section, a very good agreement with the experimental σ 0 ( ) line is observed for BW3 and BW2 models, while a slightly poorer result is provided by the PAN model.In the inset of Figure 8b, the results of the Voronoi FE model are furthermore compared with the results of the RVE model that was used in the calibration procedure.Practically identical responses observed for all three hardening models prove that the strategy employed in Section 3.2 for the identification of the RVE model was appropriate.In Figure 9, the Voronoi model was used to compare the local stress fields, in terms of von Mises stress σ mis , among the hardening models.The comparison was performed at nominal strain 0.2, where average tensile stress is roughly σ 33 = 500 MPa.As shown in Figure 9, local stresses are largely driven by the geometry of the grains, which develop stress concentrations close to grain boundaries with amplitudes up to ∼2 times the average stress.In such locations, also the influence of different hardening models is expected to be more pronounced due to amplification.However, local differences between BW3 and BW2 models are practically invisible in Figure 9, but become considerably larger when compared to the PAN model (note ∼40 MPa larger average stress in PAN).With PAN model, for example, relatively higher stress concentrations are predicted at several grain boundaries, as shown in Figure 9.  8a) at nominal strain e = 0.2 using the three hardening models with best-fit parameters P 2 , P 6 and P 10 from Table 4.
In order to quantify the differences between the BW3, BW2 and PAN models the subtracted von Mises stress fields from Figure 9 are presented in Figure 10a.As anticipated, some of the biggest variations of von Mises stress ∆σ mis are observed close to grain boundaries, however, with no clear correlation between ∆σ mis and σ mis .An obvious difference of the involved stress ranges is visible in Figure 10b: the standard deviation of ∆σ mis between the BW3 and BW2 models is estimated as ∼5 MPa (i.e., 1% of σ 33 ), while an order of magnitude bigger value of ∼70 MPa (i.e., 13% of σ 33 ) is observed between the BW3 and PAN models.Although macroscopic differences between the latter two models are relatively small (e.g., compare the fits in Figures 5b and 7b), substantially larger local stress deviations imply that great caution should be taken when interpreting the results of local fields.It seems, in fact, that grain-scale CP model calibration against tensile curves may not be sufficient to validate the model also on a local (sub-grain) scale.This is especially relevant, for example, when performing polycrystalline FE simulations of (intergranular) crack propagation (e.g., [6,19,21]), where strong geometrical amplification of stresses close to a crack tip is expected to increase the sensitivity of the model response (i.e., crack growth) to the involved uncertainties in the CP model.

Conclusions
The purpose of this work was to demonstrate a strategy for reliable calibration of the conventional, length-scale independent crystal plasticity model when using combined experimental data from the meso (single crystalline) and macro (polycrystalline) scales.In particular, the need for proper adjustment of the polycrystalline tensile data was emphasized and demonstrated by subtracting the length scale effect, due to grain boundary hardening, according to the Hall-Petch relation.The introduced automated calibration procedure was shown to be effective when the computationally-expensive polycrystalline finite element model was replaced by the small, but representative volume element model of the microstructure.The procedure for its identification was presented and validated by comparing with larger Voronoi-like aggregates.Finally, a simple hardening model upgrade was proposed to incorporate the grain size effects in conventional crystal plasticity.
The constitutive model was calibrated and validated with the published experimental data for 316L stainless steel obtained at room temperature.The model was found to be very accurate in reproducing experimental results, implying that the Bassani and Wu hardening law is accurate enough to predict tensile behavior of 316L stainless steel at both single-and poly-crystalline scales.However, strong geometrical amplification of local stresses close to crack tips or grain boundaries is expected to increase the sensitivity of the behavior response to involved uncertainties in the crystal plasticity model; therefore, great caution should be taken when interpreting the results of local (sub-grain) fields.
In addition, the following main findings about the deformation behavior of 316L stainless steel at room temperature were gained in this work: • A relatively small value of γ I 0 ∼ 0.02 and a high value of f 0 ∼ 0.67 in the Bassani and Wu hardening law indicate a short Stage I region and hard activation of cross slip during Stage II hardening.A low value of q ∼ 0.2 also suggests that self-hardening dominates over latent hardening.
• High similarities between two-stage and three-stage Bassani

Figure 1 .
Figure 1.(a) Typical τ α − γ α curve of a pure FCC single crystal loaded in a uniaxial tension with an initial orientation for single slip.(b) Polycrystalline tensile data of 316L stainless steel at room temperature for various average grain sizes D .Data were extracted from [31].The red line denotes the extrapolated tensile curve σ 0 ( ) using the Hall-Petch relation.

8 Figure 2 .
Figure 2. Examples of polycrystalline aggregate models with realistic gauge geometry: (a) N g = 64, N eg = 8, (b) N g = 64, N eg = 64, and (c) N g = 512, N eg = 64.(d) Single crystal model with N eg = 48.The colors denote different grains with common crystallographic orientation.Grains are meshed by linear hexahedral elements (C3D8) in polycrystalline models and by quadratic hexahedral elements with reduced integration (C3D20R) in a single crystal model.Arrows in (a,d) denote the boundary conditions used in tensile simulations.

33 *Figure 4 .
Figure 4. Convergence analysis of tensile stress σ 33 calculated at (a,b) 33 = 0.1 and (c,d) 33 = 0.2 as a function of the number of grains N g and the number of elements per grain N eg .Dashed lines represent a fit with a power law; Equation (9).Extrapolated stress value σ * 33 is marked by a horizontal dashed line.The representative volume element (RVE) model is marked by a red square.

Figure 5 .
Figure 5. Results of the two-scale calibration showing a comparison between calculated tensile curves of the crystal plasticity (CP) model with the Bassani and Wu hardening law (BW3) hardening model and corresponding measurements on 316L stainless steel [30,31].The calculated lines are fitted to the experimental data using four different fitting domains: (a) e max = 0.1, (b) e max = 0.2, (c) e max = 0.3 and (d) e max = 0.4.The corresponding hardening parameters are shown in Table4.

Figure 6 .
Figure 6.(a) Results of the calibration of the modified BW3 hardening model, Equations (14) and (15), calculated with newly identified parameters k 0 and k 1 (shown in legend in units of MPa √ µm)

Figure 7 .
Figure 7. Results of the two-scale calibration showing a comparison between calculated tensile curves and corresponding measurements on 316L stainless steel [30,31].The calculated lines are fitted to the experimental data on a fixed strain domain e max = 0.2 using four different strategies: (a) BW2 hardening model, (b) PAN hardening model, (c) BW3 hardening model, but with raw polycrystalline data with D = 3.1 µm, and (d) BW3 hardening model, but with raw polycrystalline data with D = 86.7 µm.The corresponding hardening parameters are shown in Tables4 and 5.

Figure 8 .
Figure 8.(a) Voronoi aggregate model with 1000 grains denoted by different colors.Grains are meshed by quadratic tetrahedral elements C3D10.(b) Tensile curves calculated with a Voronoi model using best-fit parameters P 2 , P 6 and P 10 from Table 4.In the inset, also the results of the RVE model are shown for comparison.

Figure 9 .
Figure 9.Comparison of von Mises stress fields calculated on a top free surface of the Voronoi model (marked by white rectangle in Figure8a) at nominal strain e = 0.2 using the three hardening models with best-fit parameters P 2 , P 6 and P 10 from Table4.

Figure 10 .
Figure 10.(a) Difference of the von Mises stress fields ∆σ mis between the BW3 and BW2 models and between the BW3 and PAN models from Figure 9.(b) Distributions of ∆σ mis calculated from all integration points of the Voronoi model at e = 0.2.

Table 1 .
A list of aggregate models used in this study.N g stands for the number of grains, N eg for the number of elements per grain and N e for the number of total elements.

Table 3 .
A list of fitting parameters with corresponding minimum and maximum values used in the calibration procedure.

Table 4 .
[40]ening parameters for 316L stainless steel obtained from the two-scale calibration procedure.Different hardening models were considered in the calibration: BW3 (Bassani and Wu[29]with Stages I, II, III), BW2 (Bassani and Wu[28]with Stages I, II) and PAN (Peirce, Asaro and Needleman[40]with Stage I).Parameters h i , τ i and χ 2 are shown in units of MPa. .

Table 5 .
[29]ening parameters for the BW3 model[29]for 316L stainless steel obtained from the two-scale calibration procedure when using raw polycrystalline data with different D and fixed e max = 0.2.Parameters h i , τ i and χ 2 are shown in units of MPa and D in µm.

Table 6 .
Hardening parameters for 316L stainless steel taken from the literature.Note that only (raw) polycrystalline data were used in the identification procedure.Parameters h i and τ i are shown in units of MPa.
and Wu hardening models indicate a very slow Stage II to Stage III transition with slip.• Grain boundary strengthening effects become relatively weak when average grain size D 90 µm, thus providing negligible additional hardening in the tensile response.