Towards the Crystal Plasticity Based Modeling of TRIP-Steels—From Material Point to Structural Simulations

With the complex multi-scale behavior of high-alloyed TRIP steels in mind, this contribution aims to complement recently established continuum mechanical modeling approaches for such materials, by considering their anisotropic inelastic response at the single crystal level. This approach generally enables the consideration of initial textures and their deformation-induced evolutions. It also represents the key theoretical and algorithmic foundation for future extensions to include phase transformation and twinning effects. Several rate-independent and rate-dependent formulations are investigated. The former are naturally associated with Karush-Kuhn-Tucker type inequality constraints in the sense of a multi-surface plasticity problem, whereas in the latter, these constraints are handled by penalty-type approaches. More specifically, the primary octahedral slip systems of face-centered cubic crystal symmetry are explicitly taken into account in our model application of the general framework and hardening models of increasing complexity are incorporated. To test the efficiency and robustness of the different formulations, material point simulations are carried out under proportional and non-proportional deformation histories. A rate-independent augmented Lagrangian formulation is identified as most suitable in the considered context and its finite element implementation as a User-defined MATerial subroutine (UMAT) is consequently studied in depth. To this end, the loading orientation dependence of the deformation and localization behaviors are analyzed through simulation of a mildly notched tensile specimen as a representative inhomogeneous boundary value problem.


Introduction
During the past decades, new steel grades with improved mechanical properties, such as high strength and pronounced ductility, have been developed, mainly motivated by light-weight applications in the automotive industry, cf. [1][2][3]. The initially fullyaustenitic steel, X3CrMnNi16-6-6, developed in the DFG Collaborative Research Center 799, clearly belongs to this group of advanced high-strength steels. Extensive mechanical and microstructural characterization, see [4][5][6][7][8][9] revealed that the mechanical properties can be attributed to the evolution of the microstructure during deformation, i.e. depending on temperature, stacking-fault energy and strain-rate, the face-centered cubic (fcc) austenite (γ ) forms twins or stacking-faults with hexagonal close-packed structure (ε) or transforms to body-centered cubic (bcc) martensite (α ). In particular, the different deformation mechanisms can occur concurrently, as shown in Fig. 24.1, which renders the formulation of constitutive models at the macroscopic scale a challenging task, cf. [10][11][12]. Although such models are already quite complex, they rarely incorporate the effect of evolving anisotropy due to texture evolution, which is of great importance in forming simulations. Furthermore, the application of such models in structural simulations is naturally associated with a length scale, at which the characteristic sizes of the microstructure are small compared to other dimensions of the problem. Therefore, employing such models to predict the behavior of devices at the micrometer scale seems to be questionable. In contrast, crystal plasticity based modeling approaches can in principle account for these effects, however, an appropriate scale-transition law has to be incorporated to give reasonable predictions at macroscopic scale. Keeping in mind that the deformation behavior of the TRIP-steel under consideration is mainly influenced by interaction of the deformation mechanisms at multiple scales-ranging from interactions between grains to interactions of stacking-faults with martensitic inclusions within a single grain-a crystal plasticity based multi-scale modeling approach seems to be even more appropriate. Although the kinematic aspects of the different deformation mechanisms are reasonable well understood, their incorporation into conventional crystal plasticity models is a challenging subject of ongoing research, especially for TRIP-steels.
Aiming for a comprehensive description of the transformation behavior in lowalloyed TRIP-steels, a material model that incorporates the stress-assisted transformation from fcc austenite to body-centered tetragonal (bct) martensite under thermomechanical loading is proposed in [13], which also takes into account the twinned martensite microstructure. The influence of the initial crystal orientation on the mechanical and the transformation behavior under homogeneous deformation is studied in [14] for the two cases of a single austenitic grain and an austenitic grain embedded in a ferritic, elastic-plastic matrix, which is described by a phenomenological, isotropic J 2 -plasticity model. In [15] the material model is extended to account for anisotropic plastic slip by means of a crystal plasticity model, which is employed both within the austenitic grains and in the ferritic matrix. A significant influence of the initial orientation of the austenitic grains-embedded in a single or oligocrystalline ferritic matrix-is found for macroscopically homogeneous deformation states.  [6] Furthermore, in [16] this model is applied to the simulation of a representative volume element of a low-alloyed TRIP-steel with an idealized microstructure containing multiple austenite grains embedded in discretely resolved polycrystalline ferrite matrix under macroscopically homogeneous thermomechanical loading. Here the influence of the sequence of thermal and mechanical loading on the mechanical response and transformation behavior is investigated. The multi-scale character of the deformation and transformation behavior of low-alloyed TRIP-steels is accounted for in [17] by assuming an idealized lath martensite microstructure and combining the elasticplastic material response of the two phases, austenite and martensite, with a criterion for the stress-assisted transformation and explicitly enforcing the compatibility and the stress equilibrium at their interfaces. In contrast to the models mentioned above, a single crystal material model for an initially fully-austenitic steel that shows a straininduced transformation from fcc austenite to bcc martensite is considered in [18] and the transformation kinetics of Stringfellow [19] is applied at the single crystal scale and the evolution of deformation bands/bands of stacking-faults that act as nucleation sites for martensite is explicitly taken into account. Good agreement between numerical simulations and experimental results from uniaxial tensile tests for a wide range of strain rates is observed. A gradient extended crystal plasticity model that includes both stress-assisted and strain-induced austenite to martensite phase transformation is proposed in [20,21] and is employed to study size effects in nanoindentation and in three-point bending tests and to investigate the influence of grain boundaries and twins in austenite grain on the transformation behavior. Consistent with thermodynamical considerations in [13,14], a single crystal model that accounts for the stress-assisted austenite to martensite transformation is discussed in [22] and extended to include twinning [23]. Here again, the influence of the initial crystal orientation on the mechanical behavior under homogeneous deformations is considered. Increasing research activity in the field of high-manganese TWIP steels has led to the development of two single crystal material models [24,25] that account for three different deformation mechanisms, namely slip, twinning and stacking-fault formation/ε-martensite formation. The latter model includes a dislocation density based hardening law, which is successfully calibrated based on quasi-static tensile tests of a polycrystaline TWIP-steel together with the corresponding microstructure evolution in terms of ε-martensite and twin volume [26]. Rather recently, a single crystal plasticity model that incorporates the stress-assisted austenite to martensite transformation is applied to the prediction of polycrystalline response of two and three phase low-alloyed TRIP-steels and the corresponding forming limit diagrams [27]. Furthermore, the well-known transformation kinetics model of Olson-Cohen [28] is extended in [29] to account for the crystallographic nature of the formation of deformation bands/bands of stacking-faults. The kinetics law is then coupled to a crystal plasticity model to describe the strain-induced transformation from austenite to martensite. The comparison between numerical simulations and the results of polycrystalline experiments shows that the temperature and stress-state dependency of the mechanical response and the transformation behavior is well captured. For further information on the application of crystal plasticity models, the reader is referred to the comprehensive reviews [30][31][32].
Although the above literature review shows that numerous models have been proposed, which account for the coupling of two or three deformation mechanisms or the different nature of stress-assisted and strain-induced martensitic transformation, a material model that thoroughly captures all the deformation mechanisms mentioned in Fig. 24.1 and their transition over a wide range of temperatures still seems to be missing. In order to develop such a model-possessing a modular structure and being robust at the same time-a numerical framework has to be identified that allows for a robust and efficient implementation. While the above mentioned models are almost exclusively based on rate-dependent formulations, it is argued in [33] that such an approach may introduce an artificial rate-dependence into the model's response, if such a formulation is chosen only for numerical convenience and not due to experimental results. Hence, a systematic study regarding different numerical implementations has to be carried out to assess the efficiency and the robustness of the corresponding stress-update algorithms also considering different model complexity. The contribution at hand aims for such a study and will provide recommendations for appropriate algorithms both for rate-dependent and rate-independent formulations. Therefore, the choice of rate-dependence in the constitutive description of the different deformation mechanisms can solely be made based on experimental observations rather than numerics. The current study employs a single crystal plasticity model that is able to capture the deformation behavior of stable austenitic stainless steels as illustrated in [34]. Thus, an adequate description of the TRIP-steel X3CrMnNi16-6-6 at temperatures T > 220 • C is considered, where dislocation glide is the main deformation mechanism, see Fig. 24.1. This model will also form the basis for further model developments, eventually providing a modular constitutive description of the TRIP-steel under consideration. In order to enlighten the effect of rate-dependence on the constitutive response and the robustness of the corresponding stress update algorithms, the study comprises four different formulations and hardening laws of different complexity.
This contribution is structured as follows. The single crystal plasticity model is discussed in Sect. 24.2, while the comparison of the different formulations under homogeneous proportional and non-proportional loading histories is considered in Sect. 24.3. In Sect. 24.4, a rate-independent formulation is employed to study the orientation dependence of the deformation and localization behavior of a mildly notched single crystal tensile specimen. Section 24.5 summarizes the main findings and outlines current and future research efforts on this topic.

Material Model
The constitutive behavior of face-centered cubic single crystals at finite deformations is described by a material model that builds on the approach proposed by Schmidt-Baldassari [35] in the rate-independent case and the formulation employed in [36][37][38] for the rate-dependent case. It is based on the multiplicative split of the deformation gradient into an elastic and a plastic part according to the proposal of Kröner [39] and Lee [40] F = F e · F p . (24.1) In addition to this kinematic assumption, the intermediate configuration defined by such a split is taken as isoclinic as suggested in [41]. The elastic behavior, defined with respect to that intermediate configuration, is assumed to be governed by the isotropic, volumetric-isochorically decoupled free energy function of compressible Neo-Hooke type The assumption of elastic isotropy is a reasonable approximation only for a few cubic metals, e.g. α-Tungsten, aluminum and vanadium [43], which are characterized by a Zener anisotropy index close to one [44]. This assumption is considered acceptable for the current contribution due to its simplicity. However, an extension of the material model to incorporate elastic anisotropy can be carried out by replacing the free energy density in (24.2) by the quadratic anisotropic free energy function discussed in [45], which is more suitable for moderately large elastic deformations. In contrast to the elastic behavior, the inelastic deformation of a single crystal is inherently anisotropic because it is governed by a finite number of distinct slip systems associated with the crystal lattice. For most face-centered cubic crystals it is reasonable to consider only the primary octahedral slip systems, consisting of {111} slip planes and 110 slip directions, see Table 24.1. The inelastic slip on the different slip systems, γ α , is linked to the plastic part of the deformation gradient and the plastic velocity gradient L p via the evolution equation in which s α and n α denote the slip direction and the slip plane normal of the system α. In addition, these vectors are of unit length |s α | = |n α | = 1 and are mutually orthogonal s α · n α = 0, where the latter property results in the inelastic part of the deformation gradient being isochoric. The onset of inelastic deformation on each slip system is described by limit surfaces of the form Herein, τ α denotes the resolved shear stress on the slip system α and is computed as while Y 0 and Y α respectively correspond to the initial yield stress and the driving force thermodynamically conjugate to the hardening variable ε α . The driving force is consequently defined as where an additional split of the free energy due to elastic and inelastic effects has been assumed. In the current contribution, two different hardening functions are considered.  Firstly, a purely phenomenological, Taylor-type hardening formulation of the form is introduced, which incorporates the cumulative inelastic slip A = α ε α and is parameterized by the asymptotic increase in the yield stress Y as well as the dimensionless shape parameter h. It can be deduced from the inelastic part of the free energy function, i.e. p,Taylor = Y A + 1 h exp (−hA) , (24.12) by means of (24.10). The application of the cumulative inelastic slip A in the hardening function (24.11) idealizes the interaction between different slip systems, but due to its simplicity, it has been extensively used to develop robust algorithmic frameworks for single crystal plasticity models [47][48][49]. Secondly, the alternative anisotropic hardening function proposed in [50], is considered, which introduces the symmetric interaction matrix h αβ in a phenomenological manner. It allows for a more complex interaction of different slip systems and contains up to 6 material constants [46]. In the current contribution the structure of the interaction matrix is adopted from [51]. The energy corresponding to this type of hardening function is formulated as a quadratic form (24.14) in terms of the auxiliary variable s α . The exponential type hardening is obtained from a non-associated evolution equation for s α , which eventually leads to The driving force Y α , introduced in (24.13), is obtained by taking the derivative of (24.14) with respect to s α and the subsequent substitution of (24.15) to eliminate s α from the resulting expression. Note that the hardening function (24.13) is exactly the one employed in the "GC model" in [36], which emanates from the small strain formulation presented in [50]. In order to close the system of equations evolution laws for the internal state variables have to be defined.
In the rate-independent formulation one obtainṡ for the hardening variables and for the inelastic velocity gradient from an associated formulation, which introduces the Lagrange multipliersλ α , that are subject to the Karush-Kuhn-Tucker (KKT) conditions Considering the limit surface (24.8) and comparing the inelastic velocity gradient given in (24.17) with the evolution equation for the plastic part of the deformation gradient (24.7), one can readily identifẏ (24.20) in the rate-independent case. In contrast to the commonly adopted, computational expensive active-set search algorithms to handle the inequalities in the KKT conditions, two different formulations of the rate-independent problem are considered here, which employ equality constraints only. Firstly, the augmented Lagrangian formulation, initially proposed in connection with crystal plasticity in [35], is employed, which takes the principle of maximum plastic dissipation as starting point and reformulates the inequality constrained optimization problem into an equality constrained optimization problem by means of so-called slack variables [52, pp. 72, 158-164]. The Lagrange multipliers are then obtained from λ α = max 0, η * α , (24.21) in which the viscosity-like parameter η * is introduced for purely numerical reasons, as it regularizes the problem, while the constraints are exactly enforced by means of the Lagrange multipliers. Secondly, a formulation based on nonlinear complementary functions (NCPfunctions) is considered. These functions have originally been proposed for constrained optimization problems, cf. [53]. A slightly more general form of NCPfunctions is introduced by Kanzow and Kleinmichel in [54] and employed in the current contribution, which reads as In particular, the parameter = 0 is chosen here and the Lagrange multipliers are obtained from (24.22) rather than from (24.18).
In the rate-dependent formulation, the KKT conditions in (24.18) are no longer applicable and the corresponding Lagrange multipliers are replaced by a potentially stress-dependent viscosity law v. This yieldṡ and Assuming that an equation analogous to (24.7) holds for the viscoplastic part of the deformation gradient F vp in the context of a rate-dependent formulation, one may identify the general format for the evolution equations aṡ To study the influence of the type of viscosity law on the deformation behavior of single crystals -in particular in the rate-independent limit -this contribution considers two specific cases: Firstly, the approach initially proposed by Perzyna [55] and extensively used in the groups of Cailletaud and Forest, cf. [36-38, 50, 51, 56-60], is considered, in which the viscosity law takes the form The positive material parameters K, n and η denote a stress-like scaling factor, the rate sensitivity exponent and a time-like parameter, respectively, where the inverse of the latter can be interpreted as a reference strain rate. The Macaulay brackets are defined as (24.28) and are identical to the max-function, m(x) = max(0, x), and the ramp function, r (x) = 1 2 (x + |x|). Secondly, the viscosity law introduced by Cuitiño and Ortiz in [61] and employed for instance in Miehe's group, cf. [47][48][49], which specifically reads is also employed in the studies presented here. Note that the current slip resistance of the particular slip system α is denoted by τ y α and possesses the same functional dependency on the hardening variable ε β as the quantity Y 0 + Y α (ε β ), but in (24.29) only α is explicitly dependent on both τ α and Y α . Therefore, the slip resistance τ y α is treated as history dependent normalization quantity, rather than an additional function of Y α . This allows one to automatically guarantee thermodynamic consistency and enforce that the slip γ α evolves identically to the hardening variables ε α .
In order to carry out material point calculations and structural simulations, the different formulations of the rate-dependent and the rate-independent material models have been implemented into the scientific computing environment MATLAB and subsequently into the finite element program ABAQUS via the User-defined MATerial interface (UMAT). The evolution equations for the internal state variables are integrated by means of an implicit Euler backward scheme and a projection technique is employed to enforce the incompressibility constraint for the inelastic part of the deformation gradient. Details of the corresponding algorithms and the associated tangent operator can be found in [62].

Material Response Under Homogeneous Deformation
In this chapter, the material model described in Sect. 24.2 is employed in the simulation of a fully deformation-controlled simple shear test as well as a non-proportional tension/compression-shear cycle. The results obtained with the different formulations are illustrated in the subsequent sections.

Simple Shear Loading
Due to its simplicity, the fully deformation-controlled simple shear test is extensively used in the literature to assess the robustness of various stress-update algorithms for single crystal plastictiy, cf. [35,47,[63][64][65][66][67][68][69]. For comparison, this test is also employed here, where the coefficients of the associated deformation gradient are prescribed as In particular, the simple shear motion, introduced in [63] is considered, in which the crystal lattice is misaligned with respect to the global coordinate axes and is characterized by the following orientation in terms of Euler angles {ϕ 1 , , ϕ 2 } = {0 • , 18.4349 • , 0 • } in Bunge notation [70], i.e. a sequential rotation about the z-, xand z-axes is considered. The material parameters employed in the simulation correspond to an ideal plastic behavior at the scale of the slip system and are specified in Table 24.2. Therefore, any effective hardening or softening observed in the subsequently shown stress-strain diagrams is attributed to the reorientation of the crystal lattice, which approaches a stable orientation for large shear γ , asymptotically leading to a constant shear stress. Additionally, the rate-dependent material response is characterized by the rate exponent n = 20 for the OM-viscosity function (24.29), while for the PCF-viscosity function (24.27) parameters are set to n = 10 and K = 10 −3 GPa, consistent with the experimentally observed range of rate-sensitivity exponents [71]. The influence of the chosen increment size γ on the stress-strain curve is depicted in Fig. 24.2. In the rate-independent formulation, the augmented Lagrangian algorithm takes 15, 60 and 600 steps to reach the final shear amplitude γ = 6, while the Kanzow NCP-function respectively requires 100, 600, 1200 and 1800 steps. It can be seen that both formulations converge to the same material response for the smallest increment size. The augmented Lagrangian formulation reproduces the characteristic features of the stress-strain curve even for very large shear increments and the stress response converges monotonically as the shear increment is refined. In contrast, the formulation based on the Kanzow NCP-function requires significantly smaller shear increments to ensure the numerical convergence of the Newton-type algorithm. But even though the numerical convergence is achieved, the simulated stress-strain curve is very sensitive to the chosen increment size, in particular in the strain range 2 ≤ γ ≤ 6. Herein, the constitutive response converges to the results of the augmented Lagrangian formulation with decreasing shear increment size in a non-monotonic manner. Although not shown in Fig. 24.2a, the same observations  have been made for the original Fischer-Burmeister NCP-function, which is recovered by the choice = 2 in (24.22).
In the rate-dependent formulation, the stress-strain curve under quasi-static loading condition,γ = 10 −5 s −1 , is determined in 60, 180 and 600 steps for the OM-viscosity function and in 420, 600 and 1800 steps for the PCF-viscosity function. Figure 24.2b shows that apart from small deviations, which are due to the different viscosity functions and viscous parameters employed, the stress-strain curves feature the same behavior up to γ ≈ 2.3. For larger shear strains however, a significant differ-ence in the computed stress response is observed. While the OM-formulation leads to a stress-strain curve possessing the same characteristics as the rate-independent case, the PCF-formulation yields a stress response for which the transition to higher stress levels is significantly delayed. The working hypothesis explaining this observation is that small differences in the evolution of the plastic slip γ α determined from the PCF-and OM-viscosity functions accumulate and that such an accumulation eventually leads to sudden deviation in the stress response at γ ≈ 2.3. The presence of small differences between the different rate-dependent formulations is clearly visible from the comparison of the normalized viscosity functions (24.27) and (24.29), as illustrated in Fig. 24.3. In order to check the hypothesis, a sensitiv-  Table 24.1. Note that, in c colored and gray solid lines represent the response of the augmented Lagrangian formulation for a perfectly oriented and a misaligned crystal, while the dotted lines correspond to the results of the rate-dependent PCF-formulation for a perfectly oriented crystal. Additionally, small black circles indicate the orientation at γ = 2.0, while the triangles denote a stable crystallographic orientation ity analysis is performed, in which the influence of a small initial misorientation (2 • ) on the stress-strain curve is investigated. The simulations are carried out by means of the augmented Lagrangian formulation of the rate-independent material model and the shear component γ is incremented in 600 equidistant steps. The corresponding results are shown in Fig. 24.4a. It can be seen that the initial part of the stress-strain curve is rather nonsensitive to the perturbation of the initial orientation and that the initial misorientation strongly affects the stress-strain curve only in the range 2 ≤ γ ≤ 4.5, which is consistent with the scatter observed in Fig. 24.2 for the different step sizes and formulations. Observing that small differences in the initial conditions lead to considerable deviations in the stress response indicates a stability problem, in particular a bifurcation. However, it is not caused by the Taylor ambiguity problem [72], where several combinations of plastic slip variables exist that lead to the same stress state. This is because, as shown in the literature, both the rate-independent augmented Lagrangian formulation as well as the ratedependent formulations yield unique solutions for the slip system selection and the corresponding rates [67,72] and thus the Taylor ambiguity problem is avoided. Conducting a systematic variation of the initial orientation according to φ 1 = φ 2 = 0 • and ∈ [0 • , 45 • ] reveals that the stability issue mentioned above is related to the discrete nature of the plastic flow of the single crystal. In fact choosing initial orientations in the range ∈ [22.5 • , 45 • ] asymptotically leads the reorientation of the crystal towards the stable orientation indicated by triangles in Fig. 24.4c. In the narrow range ∈ [17 • , 22 • ] the crystal approaches a different stable orientation, while for ∈ [0 • , 16 • ] a constant rotation of the crystal lattice is predicted by the material model. As the initial orientation {ϕ 1 , , ϕ 2 } = {0 • , 18.4349 • , 0 • } and the perturbed initial orientations considered are very close to the boundary of the two ranges of , the sensitivity of the stress response to small perturbations is not surprising. It is therefore concluded that small initial misorientations or the accumulation of small differences in the evolution of the plastic slip variables in the rate-dependent case are responsible for selection of slip systems made by the corresponding algorithms. In addition to the stress-strain response, Fig. 24.4b illustrates the evolution of the hardening variables for three different cases. In particular, a rate-independent, perfectly oriented crystal as well as a rate-independent crystal with a specific misorientation of 2 • , indicated by an arrow in Fig. 24.4a and a rate-dependent, perfectly oriented crystal based on the PCF-formulation is considered. Note that the slip systems, which are active only in the range 0 ≤ γ ≤ 1 and overall show a low activity, i.e. ε α < 0.35, are not shown for the sake of clarity. Comparing the pattern of activation and deactivation of the slip systems of the perfectly oriented and the misaligned crystal, it is obvious that the initial misorientation delays the deactivation of the slip systems D4 and D6 as well as the activation of the systems C3, C5, A2 and B2 (15) , where the latter corresponds to slip system 15 in Table 24.1. This delay is responsible for the postponed transition to the orientation {ϕ 1 , , ϕ 2 } = {90 • , 45 • , 0 • }, for which a constant yield stress and no rotation of the crystal is observed [63]. Furthermore, it is observed that the pattern of activation and deactivation of the slip systems computed for the perfectly oriented, rate-dependent crystal is similar to the pattern obtained for the misaligned, rate-independent crystal. This indicates that the slip γ α deter-mined via the rate-dependent formulation causes a small misalignment of the crystal compared to the results of the rate-independent response of the perfectly oriented crystal employing the augmented Lagrangian formulation, as shown in Fig. 24.4c. This leads to a different reorientation of the 100 directions and a delayed transition towards the orientation {ϕ 1 , , ϕ 2 } = {90 • , 45 • , 0 • }. These observations therefore support the hypothesis introduced above and give an explanation for the differences in the stress response presented in Fig. 24.2.
Technically, the misorientation R is incorporated by means of the Euler-Rodrigues formula (24.31) where is the Levi-Civita tensor, while a denotes the rotation axis and α the corresponding angle. The misorientations employed in this section are generated by taking the rotation axis according to the 42 equally spaced points on the unit sphere [73] and the rotation angle α = 2 • . The initial orientation of the crystal with a predefined misorientation R * is then obtained, according to [74, p. 68] from in which R is the unperturbed initial orientation.

Non-proportional Tension/compression-Shear Loading
While the fully deformation-controlled simple shear test provided valuable insight into the robustness of the different formulations in the ideal plastic case, the capabilities of the different formulations in the hardening case are assessed now under complex loading conditions employing a non-proportional load cycle, similar to the one proposed in [38]. It corresponds to a uniaxial tension/compression loading combined with a simple shear loading and approximately mimics the deformation path observed in the dual actuator loading system described in [75] or a thin-walled tubular specimen in a tension-torsion testing device [76]. The temporal change of the coefficients of the deformation gradient tensor associated with this so-called butterfly test are prescribed according to and are illustrated in Fig. 24.5a. The remaining coefficients, indicated by * , are determined by an iterative procedure, the constitutive driver, described in [77,78] that enforces the associated components of the first Piola-Kirchhoff stress to vanish, see also [79, p. 562f] for a spatial formulation. Furthermore, this procedure is enhanced by an adaptive time-stepping algorithm to allow for a step size adjustment of the prescribed components F 12 and F 22 based on the convergence behavior.
Recall that the degree of anisotropy in the hardening behavior of the single crystal material model depends on the choice of the interaction matrix in the hardening law, given in (24.13). Since its coefficients can be estimated from non-proportional tests involving strain path changes, cf. [80], it is interesting to study the influence of the interaction matrix on the performance of the different formulations in case of the butterfly test. Therefore, the initial orientation of the crystal is taken as {ϕ 1 , , ϕ 2 } = {0 • , 0 • , 0 • }, i.e. the crystal axes are aligned with the axes of the global coordinate system and the hardening parameters are chosen according to Table 24.3, while the remaining parameters are taken from Table 24.2. The coefficients of the interaction matrix are adopted from the literature without modifications, while the hardening parameters Y are rescaled in the cases (ii) and (iii). This rescaling is necessary in order to solely study the effect of the deviation of the interaction matrix from the Taylor-type hardening case, because according to (24.13) keeping the parameter Y constant would result in considerably different hardening rates. Therefore, the hardening parameters Y are adjusted in such a manner that the stress-strain curves coincide under simple shear loading for the arbitrarily chosen orientation The calibrated material model is used to simulate the stress response to the non-proportional deformation path illustrated in Fig. 24.5a employing the rateindependent, augmented Lagrangian formulation and the rate-dependent formulation based on the OM-viscosity function, because both formulations proved robust and yielded consistent results in the fully deformation-controlled test conducted in Sect. 24.3.1. For the rate-independent formulation, the computed stress path is illustrated in Fig. 24.6a, in which the reference simulations, depicted with colored lines, required 1000, 2000 and 4500 increments for the different interaction matrices in cases (i), (ii) and (iii), respectively. Furthermore, white circles denote simulation results obtained with coarser increment sizes in the cases (i) and (ii), in which 50 and 1500 increments were employed, still yielding sufficiently accurate predictions of the stress path. It is clearly visible from Fig. 24.6a that the interaction matrix significantly influences the stress response, especially at the latter stages of the deformation path, suggesting that the material parameters entering the interaction matrix could in principle be identified from such a non-proportional test.
Besides this observation, the choice of interaction matrix also significantly influences the performance of the augmented Lagrangian formulation in conjunction with the constitutive driver. In particular, while a rather coarse incrementation can be chosen in case (i), i.e. an interaction matrix corresponding to Taylor-type hardening is employed, the beneficial robustness of the formulation is diminished in case of the two-parameter interaction matrix (ii) and especially in case of the four-parameter interaction matrix (iii), at least if the same accuracy as in the reference simulation is sought. This aspect is also illustrated in Fig. 24.6b, which shows the evolution of the norm of prescribed deformation gradient components for the reference simulations. While no adaptive adjustment of the increment size is necessary in cases (i) and (ii), it is worth noting that at least one and a half times the number of increments are required in case (ii) compared to the reference simulation with Taylor-type hardening (i) to obtain a stress response independent of the chosen step size. In case (iii), an adaptive adjustment of the increment size is even necessary, particularly at the change in the deformation path from the combined tension/compression-shear loading to only shear loading in order to ensure convergence of the constitutive driver, yielding increment sizes below 10 −6 . Thus, although the augmented Lagrangian formulation proved very robust in a fully deformation-controlled test in the absence of hardening, similar results cannot be obtained in situations where the hardening is anisotropic. This is due to the inclusion of interaction matrices other than the one corresponding to Taylor-type hardening and in particular in situations that involve iterative procedures to ensure stress-free conditions in certain directions. This is not surprising, because interaction matrices not associated with Taylor-type hardening allow for a change of the shape of the elastic domain and not only its size, which confronts the corresponding stress-update algorithm with a significantly more difficult task.
Furthermore, the rate-dependent formulation based on the OM-viscosity function is now employed in the simulation of the non-proportional tension/compression-shear cycle, where the deformation gradient coefficients are prescribed atḞ 22 = 10 −2 s −1 andḞ 12 = 5 × 10 −3 s −1 , respectively. The computed stress path is illustrated in Fig. 24.7a. Again, reference solutions are indicated by colored lines. Here, they were respectively obtained with 1000, 2000 and 3300 increments for the three different interaction matrices. In the cases (i) and (ii), additional simulations have been carried out with a coarser incrementation, employing 800 and 1000 increments. The corresponding results are again depicted with white circles. For these two cases, only minor differences are observed compared to the rate-independent results and the deviations can be attributed to the chosen viscosity parameters and the moderate loading rate. Note that the stress-update algorithm corresponding to the rate-dependent OMformulation allows for a further reduction of the number of increments in case (i) and (ii), but this results in a significant loss of accuracy in the predicted stress path.
In case (iii), associated with the four-parameter interaction matrix, the stress path predicted by the rate-dependent formulation differs significantly from the results of the rate-independent formulation. A noticeable difference is already observed in the first part of the loading cycle, corresponding to a constrained uniaxial tensile loading, in which F 12 = F 21 = 0. Here, non-zero shear stresses σ 12 of considerably different magnitude develop for both the rate-dependent and rate-independent formulation, as depicted in Fig. 24.8b. These stresses can be attributed to the above mentioned constraint in combination with the activation of the slip systems given in Table 24.4. While the four-parameter interaction matrix provokes the activation of five slip systems in the rate-independent formulation, eight slip systems are active in the rate-dependent formulation, as shown in Fig. 24.8a. Note that for the latter the eight hardening variables ε α do not evolve identically. On the contrary, the simulations conducted with the interaction matrices (i) and (ii), identical slip along the same set of active systems is determined. Thus, the high-symmetry of the initial orientation is preserved during the first part of the loading cycle and the constraint is automatically fulfilled.
Comparing the complete stress paths obtained for the rate-dependent and the rateindependent formulations in case of the interaction matrix (iii) once again, it becomes apparent that the different slip activity during the initial stage of the non-proportional test results in a more complex stress path for the rate-dependent formulation. This observation is also reflected in a different history of the prescribed deformation increments shown in Fig. 24.7b, which is obtained by the adaptive procedure mentioned above. However, the sudden decrease in the increment size, indicated by a black cross in Fig. 24.7a, b, is not linked to any sharp change in the stress path. It is rather caused by rapidly changing shear components of the deformation gradient, determined from corresponding zero stress condition, to allow for a constant slip activity. Similar to the approach employed in Sect. 24.3.1, a sensitivity analysis is conducted to investigate the deviation of the stress paths during the initial constrained uniaxial tensile loading. To this end, the initial orientation is perturbed employing (24.31) and (24.32) and choosing α = 0.5 • . A wide range of stress paths is obtained for selected perturbed initial orientations employing the rate-independent augmented Lagrangian formulation. The corresponding results are illustrated in Fig. 24.8b together with the uniaxial stress path of the perfectly oriented single crystal simulated employing the rate-dependent OM-viscosity function. The striking similarity in the results of the two formulations is the sudden change of the shear stress σ 12 , which is observed immediately at the initiation of plastic flow in the rate-independent formulation for some perturbed initial orientations. In the rate-dependent result however this sudden deviation from uniaxial tension is shifted to higher stresses.
Due to the fact, that the loading axis is aligned with a direction of high symmetry for the perfect cube orientation {ϕ 1 , , ϕ 2 } = {0 • , 0 • , 0 • }-yielding up to eight potentially active slip systems-the stress response obtained either experimentally or from numerical simulations is rather sensitive to small misorientations [61,68]. As mentioned above the interaction matrix in case (iii) initially provokes the activation of this set of eight slip systems associated with the high symmetry orientation in the rate-dependent formulation (see Fig. 24.8a and Table 24.4). However, at an uniaxial stretch F 22 ≈ 1.02, several slip systems are deactivated and concentrated slip on system B5 is observed, which accompanies the sudden change in the shear stress σ 12 in the response of the rate-dependent model.
The strong sensitivity of the actual stress path to small initial misorientations, illustrated in Fig. 24.8b, indicates an unstable orientation, which is also supported by two results found in the literature. Firstly, the combinatorial search conducted in [67] for a single crystal-with an anisotropic hardening law (based on a six-parameter interaction matrix)-under incompressible, uniaxial tension loading in cube orientation revealed that there exist three different slip system solutions. One corresponds to the activation of eight systems, while the other two only activate four systems. This is consistent with the results presented in [72], where it was also found that the Taylor ambiguity problem occurs if anisotropic hardening is included in the model for fcc single crystals. Secondly, it has been reported in [72] that the value of the rate-sensitivity exponent influences the stability of a crystallographic orientation. In particular, crystallographic orientations which are stable in simulations carried out by rate-independent formulations can become unstable in the corresponding ratedependent formulation, if the viscous parameters are not chosen such that they are able to recover the rate-independent limit. The choice of viscous parameters and the Taylor ambiguity problem are therefore regarded as the key factors for the significant deviations in the stress path and slip activity obtained by the rate-independent and rate-dependent formulations. The large difference in the plastic slip activity may also be amplified due to application of the iterative procedure to determine the non-zero shear components of deformation gradient, which eventually leads to an unsymmetric plastic slip. The initial variation in stress path clearly has an influence on the subsequent evolution of the stress along the remaining non-proportional deformation path. Thus, it cannot be expected that the stress path obtained by the rate-dependent formulation in Fig. 24.7 returns to the one of the augmented Lagrangian formulation in Fig. 24.6 for the interaction matrix employed in case (iii).

Constrained Tension Test
To further illustrate the robustness of the rate-independent, augmented Lagrangian formulation, it is tested within the finite element solution of a spatially inhomogeneous finite deformation boundary value problem. In particular, a constrained tension test is simulated employing a tensile specimen with a geometrical imperfection. The corresponding dimensions are given in Fig. 24.9, where the imperfection is introduced in terms of two symmetrical notches reducing the width of the specimen at its center to 90% of the original width. Along the lower edge of the specimen the displacement in 2-direction is fixed, while at the upper edge a proportionally increasing displacement is homogeneously prescribed. Additionally, the influence of two different types of boundary conditions is studied. Firstly, a clamped condition, in which the displacement in 1-direction u 1 = 0 is prescribed at both edges and secondly, a free lateral contraction condition, for which the displacement in 1-direction u 1 = 0 is only enforced at points A and B in Fig. 24.9, is considered. In either case, the displacement in 3-direction is also fixed for the entire specimen. The geometry of the specimen is discretized with 400 three-dimensional linear hexahedral elements, denoted as C3D8 according to ABAQUS conventions. The material model outlined in Sect. 24.2 with Taylor-type hardening, (24.11), is employed in the simulation. The corresponding material parameters are given in Table 24.5. With this model at hand, the influence of the initial orientation of the crystal on the force displacement curve is studied first. Therefore, four different initial orientations with respect to loading axis, namely {ϕ 1 , ,      Fig. 24.10a. For the orientation {289 • , 163.6 • , 42.5 • }, a considerably lower force is required at initial yield, but the plastic flow is accompanied by strong initial hardening, leading to almost the same force level as the former orientation ( Fig. 24.10a).
Besides the influence of the initial orientation on the initial yield and the hardening behavior, the onset of necking and localization of the deformation is also strongly affected by the choice of initial orientation. While the force-displacement curves for the orientations {0 • , 0 • , 0 • } and {0 • , 45 • , 30 • } show a rather smooth transition to the geometrically-induced softening, more rapid drops in the applied force are observed for the other two orientations. A similar trend is observed for the respective minimal width of the specimen w min measured in 1-direction, shown in Fig. 24.10b, where the former two orientations continuously deviate from the initially linear relation between w min and u 2 , while especially the orientation {30 • , 45 • , 0 • } is characterized by a sudden change in the specimen width. In each case, the deformation at which the deviation from this linear relation occurs, does not coincide with the maximum applied force, but happens, as expected, at considerably smaller deformations. The deformed specimen geometry and the spatial distribution of the logarithmic strain in longitudinal direction, i.e. (LE 22 ) are illustrated in Fig. 24.12, where the latter is computed from The evolution of the strain fields further confirms the important influence of the initial orientation on the mode of localization (symmetric/unsymmetric). In case of the highly symmetric crystal orientation, {0 • , 0 • , 0 • }, the strain field is symmetric even during necking. The location, at which the highest strain is observed, is shifted rapidly from the notch surface, for the elastic solution, shown on the left in the top row of Fig. 24.12, to the center of the specimen during elastic-plastic loading. In contrast, the strain field of the initial orientation {289 • , 163.6 • , 42.5 • } shows a smooth transition at the initiation of plastic deformation. While during the initial stages of elastic-plastic loading a rather large volume experiences pronounced deformation,  the deformation eventually localizes within a band. Interestingly, the inclination of the deformation band changes during the loading history due to the reorientation of the crystal. Motivated by similar studies presented in [81], the influence of the boundary conditions indicated in Fig. 24.9 shall now be analyzed. Inspecting the respective forcedisplacement curves, depicted in Fig. 24.11a, only minor differences are observed for these two scenarios. In particular, the clamped condition leads to a slight increase of the required force at moderate elongations (u 2 ≤ 5 mm) and to an earlier localization, compared to the free lateral contraction condition. Especially, the latter observation is also clearly visible in the evolution of the minimal specimen width, see Fig. 24.11b. Although the influence of the boundary condition is not as pronounced as reported in [81], the main characteristic, namely the earlier localization for the clamped condition is consistent with the results from the literature. The reason why the current simulation only shows a small sensitivity with respect to the boundary condition is that the reference case employed a double slip formulation, while the current model includes all primary slip systems of the fcc material. Moreover, inspecting the strain fields in Fig. 24.13 reveals that the clamped condition has significant influence only in the initial loading stages, where it leads to the development of a localization band immediately after plastic yielding-an effect that is absent in the free lateral contraction case. However, as the prescribed deformation is increased the strain fields become increasingly similar, although they differ in absolute values of LE 22 .
Finally, it is worth emphasizing two aspects that pertain to all of the numerical tensile test studies. Firstly, although significant deformation increments can occur at the Gauss point level for elements within a developing deformation band, the material routine, in fact, did not require a reduction in the global time step. This again emphasizes the robustness of the augmented Lagrangian formulation, also in the context of inhomogeneous deformation states. Secondly, it is well-known that the formation of localization bands in the geometrically-induced softening regime leads to spurious mesh-dependencies of numerical results. The reason is that in local formulations-which lack the notion of an intrinsic length scale -localization zones degenerate to discontinuities surfaces, whose predicted widths solely depend on the spatial resolution of the finite element mesh. Generalized continuum formulations have been proposed in the literature that circumvent this problem and also naturally incorporate size effects, cf. [82][83][84][85] for gradient extended formulations and [86][87][88][89] for micromorphic formulations. However, the application of the proposed material model within the framework of a generalized continuum formulation is beyond the scope of the current contribution.

Conclusions
In the current contribution, a material model for face-centered cubic single crystals, suitable for finite deformations, was discussed. Four different formulations were considered accounting for both rate-dependent and rate-independent flow behavior. Furthermore, a nonlinear anisotropic hardening law based on an interaction matrix has been incorporated, which accounts for the interaction of different slip systems in a phenomenological manner. The robustness of the corresponding stress update algorithms was assessed under homogeneous deformation states employing proportional and non-proportional loading histories. In the rate-independent case the augmented Lagrangian based formulation, originally proposed in [35], proved to be very robust, while the rate-dependent model adopted from [61] provided convincing results. However, increasing the complexity of the hardening law by choosing the parameters in the interaction matrix increasingly different from Taylor-type hardening, had a strong impact on the performance of the corresponding algorithms, resulting in a strong reduction of the prescribed deformation increment size. But this observation can readily be explained by the fact that such anisotropic interaction matrices induce an evolution of the shape of the elastic domain and not only its size, complicating the stress update considerably. Thus, the successful application of these two stress-update algorithms even in the case of anisotropic hardening emphasized their robustness.
In order to assess the performance of the rate-independent, augmented Lagrangian formulation under inhomogeneous deformations, an implementation of this model as a user-defined material subroutine (UMAT) into the finite element program ABAQUS has been employed in the simulation of a single crystalline, mildly notched tensile specimen. Herein, a strong influence of the initial orientation on the deformation and localization behavior was observed. Even during the formation of localization bands in the specimen-resulting in substantial deformation increments at the Gauss point level-an adjustment of the global time step was not necessary, confirming the robustness of the augmented Lagrangian formulation in connection with a Taylortype hardening law also under inhomogeneous deformation states.
Having identified robust numerical frameworks for both the rate-dependent and the rate-independent case, current research efforts are focused on model extensions towards the inclusion of martensitic phase transformation and twinning by means of analytical homogenization approaches. Furthermore, a comparison with experimental results at the single crystal scale is sought, from which the constitutive parameters will be identified.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.