Solving Material Mechanics and Multiphysics Problems of Metals with Complex Microstructures Using DAMASK—The Düsseldorf Advanced Material Simulation Kit

Predicting process–structure and structure–property relationships are the key tasks in materials science and engineering. Central to both research directions is the internal material structure. In the case of metallic materials for structural applications, this internal structure, the microstructure, is the collective ensemble of all equilibrium and nonequilibrium lattice imperfections. Continuum models to derive process–structure and structure–property relationships are based on two ingredients: 1) quantitative state variables that capture the essential features of the material's state and 2) kinetic equations for the state that describe its evolution under load. Successful models, that is, models that are of practical use, depend on state variables and corresponding evolution laws that are sufficiently representative for the microstructure and that are able to describe the phenomena of interest. The development of software tools capable to integrate these different aspects to get a holistic view of process–structure–property relationships requires joint efforts from specialists in different disciplines and a long‐term perspective. The Düsseldorf Advanced Material Simulation Kit (DAMASK) is such a tool. In this overview article published on the occasion of Advanced Engineering Materials' 20th anniversary, some representative application examples which demonstrate how DAMASK can be used to study metallic microstructures at different length scales are presented.


Introduction
The prediction of process-structure and structure-property relationships are the two key research fields in (computational) materials science and engineering. [1][2][3] Central to both areas of research is the internal material structure. In the case of metallic materials used for structural applications, this internal structure is called microstructure and is the collective ensemble of all equilibrium and nonequilibrium lattice features including 0D, 1D, 2D, and 3D defects. While investigations on processstructure relationships aim at answering how a thermo-chemo-mechanical process influences the microstructure, investigations on structure-property relationships aim at answering how a given microstructure responds to applied thermo-chemomechanical loads. Connecting process, structure, and properties in a systematic way requires to 1) arrive at a description of the microstructure and 2) identify how it reacts to thermal, chemical, and mechanical loads. Examples for microstructure descriptors are dislocation density, phase volume fractions, and crystallographic texture; examples for loads are controlled tensile experiments, industrial cold rolling, or sour gas environments during service.
Developing continuum models to derive process-structure and structure-property relationships can be divided into two tasks: 1) the identification of adequate quantitative state variables that represent the most essential features of the material's microstructural state, substructure, and texture 2) the establishment of kinetic equations for the state that describe the evolution of these state parameters under applied boundary and/or initial conditions. The art of successful modeling lies in selecting state variables and corresponding evolution laws that are able to describe the phenomena of interest. Examples for state variables and their respective evolution laws are the dislocation density according to the Kocks-Mecking law, [4] the chemical composition that evolves according to Fick 's law [5] or according to the Ginzburg-Landau theory [6] in spinodal decomposition, [7] or the order parameter in phase field models. [8] The level of detail required for a suitable microstructure representation is closely related to the focus of the investigation: Although a single dislocation density variable is often sufficient when investigating the collective behavior and strength of an averaged ensemble of multiple grains, investigation of grain-grain interactions might make it necessary to further separate the dislocation density according to their screw and edge character [9][10][11] or whether they are statistically stored or geometrically necessary. [12,10,13] For practical reasons, it is advisable to keep the number of the internal variables and constants of the evolution equations as small as possible to capture the relevant effects and select them such that their values are accessible to experimental or computational testing protocols. Finally, as the models are typically embedded into a discretized continuum description, the resulting set of equations should be solvable by an efficient and robust numerical method.
Despite the close connection between process-structure (such as recrystallization rate in dependence of temperature profile) and structure-property (such as strength in dependence of dislocation density and texture) relationships, both aspects are often investigated separately and with different models and numerical tools. On the one hand, investigations on process-structure are performed with phase field methods, [8] cellular automata, [14] or vertex (front-tracking) approaches. [15] Focusing on the effect of temperature and composition on moving internal interfaces, the fact that metallic materials are often subjected to heavy deformation during processing, e.g., in the case of hot rolling, is often neglected. On the other hand, structure-property relationships are investigated with the help of constitutive equations for continuum damage evolution [16] or crystal plasticity [17,18] and finite element method (FEM)-or spectral/fast Fourier transform-(FFT)-based solvers. [19] Focusing on the deformation caused by mechanical loads, one often neglects the influence of chemical or thermal loads during service that might decrease the lifetime-for instance, in the case of hydrogen embrittlement-significantly. The mutual thermo-chemomechanical interactions that determine both, the evolution of the microstructure and its response to loading during production and in service, however, motivates the development of coupled multiphysics modeling frameworks. Such frameworks enable to determine process-structure-property relationships or, from an industrial perspective, serve as an Integrated Computational Materials Engineering (ICME) tool capable to simulate the materials chain from production to service.
The development of simulation tools for microstructurecentered investigations has seen a large progress since the establishment of Advanced Engineering Materials (AEM) in 1999. Comparing the work from the early volumes [20][21][22][23][24][25][26] with recently appeared publications [27][28][29][30][31][32][33][34][35] reveals not only that the increased computational power allows to increase time and length scale and temporal and spatial discretization, but also that recent models are more often based on physical-oriented grounds instead of relying on phenomenological expressions. [36] A rather recent trend, only made possible by these improvements, is the integration of models describing different aspects such as dislocation glide, twinning, phase transformation, recrystallization, and precipitation into truly integrated tools for ICME. [37] In contrast to the often more quantitative progress, combined multiphysics tools introduce a qualitative step forward in materials modeling as they allow to investigate how different microstructural mechanisms of responding to an applied load can either compete with or mutually enhance each other. Such integrated tools allow to investigate, for example, if a material realizes the prescribed deformation by crystallographic slip or mechanical twinning, [38] how the plastic and elastic anisotropy determines the morphology of second-phase particles during precipitation, [39] whether a material hardens by deformation or softens by recrystallization, [40][41][42][43] and whether plastic deformation or fracture is energetically more favorable under a given loading scenario. [44] The development of tools capable to integrate these different aspects to get a holistic view of process-structure-property relationships requires joint efforts from specialists in different disciplines and a long history in sustainable software development.
The Düsseldorf Advanced Material Simulation Kit (DAMASK) [45,46] is such a software simulation tool for ICME. DAMASK is available as free and open-source software and can be downloaded at https://damask.mpie.de. Originally intended to be a crystal plasticity subroutine, that is, focusing more on structure-property relationships (making its first appearance in AEM in ref. [47]), it has emerged into a modular package which is designed in a hierarchically structured way of modeling the material point behavior for the solution of elastoplastic boundary value problems along with damage, thermal and some thermodynamical effects from the single crystal up to the scale of parts in complex manufacturing processes. It has been applied to a wide range of engineering alloys, focusing on crystal plasticity and multiphysics problems. Models range from phenomenological descriptions to physics-based formulations of dislocation slip, [11] twinning-induced plasticity (TWIP), martensitic transformations (TRIP), [48] microstructural damage evolution, [44,49] and the associated temperature evolution. [19] The resulting sets of nonlinear internal-variable differential equations are solved in a fully coupled way using either FEM or spectral-based FFT solvers. DAMASK is used by multiple academic and industrial users worldwide. [50][51][52][53][54][55] In this overview article published on the occasion of AEM's 20th anniversary, we present recent application examples showing how DAMASK can be used to investigate the response of metallic microstructures to applied loads. While, due to DAMASKs origin as a crystal plasticity framework, most examples presented here deal with the materials response to mechanical loads, they also demonstrate-along with several recent publications [19,[57][58][59][60] -its capabilities as a multiphysics ICME tool. The examples presented here are mainly based on existing publications and include yield-surface-based forming simulations, Lankford value predictions, and the forecast of damage in multiphase steels.

Multiphysics at and below the Grain Scale
Mechanism-oriented alloy design warrants the understanding of strengthening and weakening mechanisms and the mutual interactions among them. Continuum modeling-which is based on existing knowledge about such mechanisms obtained from experiments or lower scale simulations-is a powerful tool to investigate under which conditions the different mechanisms are activated and how they interact. If the probed volume is large enough to contain all relevant microstructure and phase constituents, that is, if the volume element is representative, macroscopic properties can be obtained from such simulations as a "side product".
The studies presented here span the scale of individual dislocations (Section 2.3 [60] ) over single grains (Section 2.1 [56] ) up to grain aggregates (Section 2.2 [38,61] ). Especially for studying the interaction of individual grains, it is a common approach to use experimentally obtained microstructures as input for the simulations. [62,63] Alternatively, tools for the creation of synthetic microstructures [64,65] are used which allow to systematically vary individual features of the microstructure. This is also the approach used to generate the microstructures for the application examples presented in the following to demonstrate how multiphysics-enhanced crystal plasticity modeling can give insights into deformation mechanisms.

Design Guidelines for Damage-Tolerant High-Modulus Steels
From a materials engineering perspective, lightweight alloy design is the development of materials with improved propertyto-weight ratio. Hence, all considered material properties should be assessed relative to the mass density ρ when evaluating the materials performance. The properties that are usually considered when rating structural materials are specific yield strength (σ y / ρ; to avoid plastic deformation during service) and specific total elongation (ε t / ρ; to enable processing and prevent detrimental failure in the case of overloading). [66,67] One important property, namely the specific stiffness E / ρwhere E is Young 's modulus-is often not considered even though an increased stiffness reduces the elastic deformation and prevents buckling. As the specific stiffness among most commercially established metallic material classes-from magnesium to steel-is virtually constant at approximately E / ρ ¼ 0.026 GPa m 3 kg À1 , increasing the stiffness-to-weight ratio imposes a severe challenge in materials design. One promising approach to design metallic materials with a high stiffnessto-weight ratio is the reinforcement of a metallic matrix with ceramic particles. The resulting material is known as a metal matrix composite (MMC). The stiff ceramic particles increase the stiffness of the MMC by up to 25% in comparison to the base material. [67] The advantage of this approach in comparison to texture design is the feasibility to increase the stiffness simultaneously in all directions. Its disadvantage, however, is the existence of numerous potential damage-initiation sites due to the presence of the brittle and stiff ceramic particles. Here, it is discussed-based on Wang et al. [56] -how coupled crystal plasticity-phase field method for fracture simulations can be used to derive guidelines for stiff yet damage-tolerant MMCs.
The material that is exemplarily investigated in this study is the iron-based Fe-TiB 2 MMC. Under suitable processing conditions, this material has a stiffness of E ¼ 240 GPa (steel: E ¼ 210 GPa) and a density of ρ ¼ 7400 kg m À2 (steel: ρ ¼ 7800 kg m À2 ). To identify favorable particle microstructures that provide optimum damage resistance, key parameters such as particle clustering degree, particle size, particle morphology, dispersion, and volume fraction [68,69] need to be systematically varied. The efforts associated with systematic experimental studies render them often impractical due to the high number of experiments and difficulties in varying multiple parameters systematically. In the studied case it is, for example, challenging to control size, morphology, and dispersion of the particles independently. The study presented here shows how these difficulties can be overcome by performing simulations on virtual microstructures with systematically varied particle parameters. This allows to derive guidelines and design goals for the further development of Fe-TiB 2 composite materials.
The 2D MMC microstructures used here were created under the assumption of uniform particle size, circular particle shape, and single crystal ferritic matrix. The particle distributions were generated by an extended random sequential adsorption (RSA) algorithm. [70][71][72] The basic principle of this approach, in more detail explained by Wang et al., [56] is the overall consideration of the particle position, the interparticle distance, the cluster size, and also the periodic boundary conditions along the geometry edges. Figure 1 shows an example of an artificially created composite microstructure. The created microstructure was discretized directly into an equidistant grid required by the FFTbased spectral solver. [19] To quantify the inhomogeneity degree of the particle distributions, the coefficient of variance of the mean near-neighbor distance (COV d ) [73,74] was derived for each microstructure where the particle neighbors were identified by Voronoi tessellations (see Figure 1c). Simulations were performed with a phenomenological crystal plasticity model coupled to a phase-field model for brittle fracture. The phase-field damage variable φ is a scalar nonlocal variable varying between 1 (undamaged state, initial situation) and 0 (fully damaged state). Details about the used phase-field model for fracture and its coupling with the crystal plasticity formluation were discussed previously by Shanthraj et al. [49] The TiB 2 particles with hexagonal crystal structure have been assumed as an elastic material with stiffness coefficients listed in Table A.1, Supporting Information. The interface energy of the TiB 2 particles, which controls the damage behavior, was estimated to be about 68 J m À2 , according to their fracture toughness and Young 's modulus. The effective damage interface thickness l 0 was taken as 0.5 μm, corresponding to twice the width of the discretization in space. A phenomenological crystal plasticity model was used for the ferrite matrix with parameters listed in Table A.2, Supporting Information. All simulations were conducted under uniaxial tension loading along the x-direction at a rate of 0.01 s À1 with a time step of 0.01 s. The spectral solver for the mechanical boundary value problem and the phase-field kinetic evolution for fracture implemented in DAMASK [49,19] were used. Figure 2 shows an example of simulation patterns on a constructed composite microstructure with heterogeneity degree COV d % 0.55. The particle fracture occurs initially at an applied uniaxial strain of 0.2461 corresponding to local strain levels of up to 1.0. Figure 1. Schematic representation of the composite microstructure in which features used to characterize the artificially created microstructures are indicated: a) particle geometry, b) microstructure of clustered distribution, and c) its Voronoi tessellation. D p : particle size; d 0 : interparticle distance; s 0 : minimum particle distance; d 1 : cluster distance; s 1 : minimum cluster distance; s 2 : standard deviation of the normal distribution within a cluster. Reproduced with permission. [56] Copyright 2018, Elsevier. 2494. In (a), the von Mises true stress σ vM distributions of the area including the earliest damage event are drawn at the loading step just before and after particle crack. The maximum and minimum stress points inside the particle before and after cracking are indicated respectively. Reproduced with permission. [56] Copyright 2018, Elsevier.
The damage area (damage variable φ ¼ 0.0) inside the particles is shown in black color. A maximum stress peak of about 6.63 GPa is found inside the particle just before the crack initiation ( Figure 2a). After cracking, the stress inside the particle decreases, especially for the damaged area where the point of minimum stress is found. This initial particle fracture event creates modulations of the local stress field, inducing adjacent load peaks and thus subsequent fracturing of neighboring particles (Figure 2b). With further loading, particle damage percolates successively to other areas (see Figure 2c, d). This particle damage process corresponds to the observation in the damage evaluation experiment of the Fe-TiB 2 MMCs. [75] For a quantitative investigation, the particle damage behavior was characterized by the initial damage strain ε D and the damage incubation strain ε Inc as functions of certain microstructure parameters. The initial damage strain ε D is the global strain value when the first point is completely fractured. The damage incubation strain ε Inc is the strain range measured from the start of the mechanically unstable stage (damage variable φ < 1.0) to the earliest complete particle fracture. Figure 3 shows the effect of particle clustering degree in terms of the variation of the initial damage strain ε D and damage incubation strain ε Inc as function of the clustering degree COV d . The homogeneous particle distributions with COV d < 0.4 show optimal performance with high values of ε D and ε Inc . However, the particle distributions with COV d > 0.7 exhibit much lower values of ε D and ε Inc , which indicates a significantly increasing preference for particle fracture in highly clustered microstructures.
Other influencing factors, such as the particle size and the particle volume fraction, were analyzed in the study which is given in full length by Wang et al. [56] Based on these investigations, the following guidelines have been derived: Fe-TiB 2 MMCs microstructures with optimum ductility and toughness should contain 7-15 vol% particles with a homogeneous distribution (COV d < 0.4). Large primary TiB 2 particles should be avoided. Figure 3. The initial damage strain ε D and damage incubation strain ε Inc as function of the particle clustering degree COV d . The increase in the clustering degree was achieved by either concentrating particles toward each other (a,b), or by shrinking the cluster size (c,d). The average interparticle spacingD 0 of the particle distribution is qualitatively represented by the height of the orange area. The reference height of the orange area, which is in the bottom of the figure, corresponds toD 0 ¼ 2D p , where D p is the particle diameter. Reproduced with permission. [56] Copyright 2018, Elsevier.

Spatially Resolved Twinning in Magnesium
Magnesium (Mg) alloys are, due to their high strength-to-weight ratio, promising candidates for the use in many engineering applications. [76] However, the wider application of wrought Mg alloys is hindered by their limited ductility at room temperature. [77,78] It is well established that the poor formability results from the contrasting critical resolved shear stresses (CRSS) between different slip and twin systems. The basal slip systems enable deformation only along the 〈a〉 directions of the hexagonal close-packed crystal structure and have typically the lowest CRSS ([18.0 AE 3.0] MPa in Mg [79] ), whereas deformation in 〈c〉 direction can be accommodated by 〈aþc〉 slip on the pyramidal systems [80,81] as well as by compression and extension twins. [82,83] As the CRSS for 〈aþc〉 nonbasal slip and compression twinning are approximately one order of magnitude higher than for basal slip systems, [84,85] for most loading conditions, only tensile twins-with a CRSS of (23.5 AE 5.0) MPa [86,79] -are frequently observed in Mg alloys at ambient temperatures. [87,88] Better understanding and corresponding modeling of the physics of tensile twinning in hexagonal materials is therefore an essential prerequisite for processing at room temperature and wider engineering use of Mg alloys. To enable such simulations, an integrated full-field dislocation density-based crystal plasticity phase-field model has been developed. [38] This model allows the concurrent description of dislocation plasticity and heterogeneous twinning behavior including twin nucleation, propagation, and growth in Mg alloys and was implemented in DAMASK. In the following, the essential aspects of the model have been compiled from Liu et al. [38] and its capabilities are demonstrated with the help of three application examples. [38,61] In short, the modeling approach consists of three elements: 1) A spatially grain-resolved crystal plasticity model is used to predict the inhomogeneous distribution of stress, strain, crystallographic texture, and dislocation density evolution; 2) A stochastic twin nucleation model to simulate f1012g tensile twins nucleation at grain boundaries in a polycrystalline microstructure; 3) A phase-field model to drive twin propagation and growth by the Landau-type relaxation of the total free energy considering both, the orientation-dependent twin interface energy as well as the elastic strain energy.
The novel combination of these three submodels enables to systematically investigate the stress redistribution and strain localization related to the twinning process and associated local deformation accommodation, such as twinning-induced slip and twin transmission across grain boundaries. Table A.4, Supporting Information, shows the set of constitutive parameters used for the examples presented in the following. Figure 4 shows results from a study by Liu et al. [38] on the influence of the grain boundary misorientation angle and dislocation slip on the twin transmission across grain boundaries in bicrystals. In the simulation, the twinning direction in the upper crystal was parallel to the applied shear direction, whereas the orientation of the bottom grain was rotated relative to the top grain to get different misorientation angles. The simulated twinning microstructure, the twinning resolved shear stress (TRSS), and the accumulated basal dislocation density at an applied shear strain of 0.6% are shown in Figure 4. It can be seen that different deformation modes are activated in the bottom grain to accommodate the shear induced by the incoming twin in the upper grain: in the case of a small grain boundary misorientation angle of 26 ( Figure 4a, top row), the twin transmits across the grain boundary into the bottom grain (outgoing twin) under the shear deformation and a TRSS concentration is observed ahead of the twin tip. In addition, a slip band along the basal plane in the bottom grain is also observed. In the case with a large grain boundary misorientation angle of 45 ( Figure 4b, bottom row), the impinged twin is blocked at the grain boundary and a distinct slip band is formed in the bottom grain. For the bicrystal with a large misorientation angle, even though the twinning-induced shear leads to a stress concentration at the intersection of the incoming twin and the grain boundary, the stress state is not capable to stimulate another twin. Furthermore, the geometric compatibility factor between the incoming twin and the basal slip in the bottom grain is relatively high (0.87), which indicates that the incoming twin-induced shear strain can be well accommodated by the activation of basal slip in the bottom grain. Therefore, the stress concentration near the intersection of the incoming twin and the grain boundary can be relaxed by the given as an inset. Adopted with permission. [38] Copyright 2018, Elsevier.
www.advancedsciencenews.com www.aem-journal.com activation of "easy" basal slip in the bottom grain, which mitigates the necessity of twin transmission. Figure 5 shows tensile twins and their interactions with dislocation slip, local stress fields, and grain boundaries in polycrystalline Mg. The simulation results reveal how tensile twins nucleate at grain boundaries or triple junctions and rapidly propagate into the grain interior. The twins impinging at grain boundaries are either blocked or transmitted across grain boundaries depending on the orientations of neighboring grains and local stress states. The TRSS distribution map shows that twin propagation leads to a severe stress redistribution ahead of the twin tips and near the transverse twin boundaries. Moreover, the inhomogeneous plastic relaxation by dislocation slip results in a heterogeneous TRSS distribution inside some of the twin bands. Even though the nominal Schmid factor for basal slip in most of the grains has relatively low values, a strongly heterogeneous distribution of basal dislocation densities is observed in both, twinned regions and parent grains. In general, a higher dislocation density is seen near the interactions of twins with grain boundaries, that is, at regions with high stress concentrations.
As a last example, the influence of precipitates on the twinning behavior-see Liu et al. [61] for further details-is studied. Precipitation hardening is a well-known mechanism for increasing the yield strength. From an alloy design perspective, it is important to realize that the individual deformation modes can be selectively triggered or blocked with different types of precipitates. Here, the influence of the precipitate size on the threshold stress for tensile twin growth is investigated. To this end, single crystals containing 4% volume fraction of spherical precipitates of uniform but varying diameters are loaded in simple shear. Figure 6a shows the twinning microstructure, dislocation densities on basal, prismatic, pyramidal slip systems, and resolved shear stress on the f1021g twin system (TRSS) during the interaction between a tensile twin and precipitates exemplary for a fixed precipitate diameter at a strain of 8.3%. It can be seen that both 〈a〉 and 〈aþc〉 slip systems are activated during the precipitate-twin interaction to accommodate the strain incompatibility between the elastically deformed precipitates and the sheared matrix. However, the plastic relaxation zone of the basal slip is much larger than that of prismatic and pyramidal slip. Figure 6a 5 shows that interaction between the twin band and the precipitates results in a very inhomogeneous stress distribution. This can be also seen from the fact that prismatic and basal slip systems are active in the vicinity of precipitate-twin intersections (see Figure 6a 3-4 ).
The stress-strain curves of the single crystal for different precipitate sizes are given in Figure 6b. Based on the results shown here, the increase in the mean threshold stress for twin growth in the presence of precipitates (Δτ CRSS ) can be derived. It is calculated as the average flow stress between the strain corresponding to the first stress peak and the final strain in comparison to the same value without precipitates (black curve). Figure 6c shows the evolution of the CRSS increment for twin growth, Δτ CRSS , and dislocation densities in dependence of the precipitate size at a shear strain of 8.3%. An increase in the CRSS increment by 54% (from 28 to 43 MPa) is seen when the precipitate size decreases from 20 to 6 elements. However, no significant change of the considered quantities can be seen for precipitate sizes exceeding 20 elements. Hence, a critical precipitate size exists. Below this critical size, the hardening can be increased by a finer dispersion of precipitates at a constant volume fraction. If the precipitate size exceeds this critical value, strengthening is only possible by increasing the total volume fraction of precipitates. The simulations also reveal that the mean back-stress in the twin band is only slightly influenced by the precipitate size (see Figure 6c). Based on these observations, Liu et al. [61] concluded that the hardening effect of smaller precipitates is governed by an Orowan-like mechanism while the back-stress mechanism is dominant for larger precipitates.
Based on the results obtained with the combined crystal plasticity phase-field model, Mg alloys with improved ductility can be designed. Microstructural parameters that can be optimized include crystallographic orientation and grain size as well as particle shape, distribution, and size. It should, however, be noted that the current version of the model is limited to one twin variant. An extension that allows to model all six tensile twin variants simultaneously is therefore planned for the future.  figure along f81 34 47g, b) resolved shear stress on the f1012g twin system, and c) basal dislocation density (logarithmic mapping) at a shear strain of 0.65%. Adopted with permission. [38] Copyright 2018, Elsevier.

Chemo-Mechanical Interactions at Dislocations in Superalloys
Given their superior high-temperature creep and fatigue behavior, Nickel-Aluminum (Ni-Al)-based superalloys are the prime choice for turbine blades in jet engines and in power plants. [89,90] Their good mechanical properties can be attributed to their single crystal microstructure, consisting of face-centered cubic (fcc) matrix (γ) and precipitate (γ') phases. The γ-matrix phase is typically a Ni-based solid solution, and the γ'-precipitate phase is chemically ordered (L1 2 ), with Al atoms at the fcc corner sites, Ni atoms at the face sites, and solute elements partitioning to these. To further improve the alloy strength, solute elements such as Wolfram (W), Rhenium (Re), Cobalt (Co), and Chromium (Cr) are often added. Although their addition is in principle aimed at alloy strengthening, solute segregation to defects such as dislocations and stacking faults [91][92][93][94] may lead to precipitate dissolution, enhanced directional coarsening ("rafting"), and degradation of mechanical properties. This application example briefly showcases results of a recently developed model for solute segregation to defects and faults in metallic alloys and its application to Ni-based superalloys. More specifically, the model is used to investigate chemo-mechanical interactions at dislocations in the ternary Ni-Al-Co system. To study (nanoscopic) defect-solute-phase interaction in multicomponent, multiphase metallic alloy systems such as Ni-Al-Co, a general phase field-based chemo-mechanical methodology [95] has recently been used to develop an atomistic phase-field chemo-mechanics (APFCM) framework which is implemented in DAMASK. In APFCM, the modeling of dislocations and faults is based on (nonconservative) order parameters associated with lattice disregistry, stacking, and chemical order. In addition, solute segregation is modeled in APFCM via (conservative) solute mass balance. Details of the atomistic phase-field model for dislocations have been previously presented by Mianroodi and Svendsen [96] and Mianroodi et al. [97] Application of APFCM to the simulation of dislocationsolute-precipitate interaction is based on numerical solution of the coupled field relations for mechanical equilibrium, dislocation motion, and solute transport in a staggered fashion with the help of finite element or spectral (e.g., FFT) methods within DAMASK. All material properties are determined in APFCM quantitatively with the help of atomistic information. In the case of the ternary Ni-Al-Co system considered here, an embedded atom interatomic potential [98] has been used to this end. The details of the APFCM model and the parameterization are reported in studies by Mianroodi et al. [60] An exemplary result is given in Figure 7, where the Co concentration in the vicinity of a γ-γ' boundary is shown. In general, the presence of solutes in the lattice results in lattice distortion and, therefore, stress. This stress drives solute segregation to dislocations and faults. www.advancedsciencenews.com www.aem-journal.com Figure 7a shows how an initially perfect 〈110〉/2 edge dislocation in γ is dissociated into leading partial and trailing partial (LP and TP) dislocations. This results in an intrinsic stacking fault (ISF) in between the partial dislocations which is decorated with Co. As the leading partial moves into the γ' precipitate, the dislocation deposits Co on the γ-γ' interface ( Figure 7b) and drags a portion of Co into the precipitate. In γ', the leading partial leaves a complex stacking fault (CSF) behind (Figure 7b). In addition, because a 〈110〉/2 dislocation does not recover chemical ordering in γ' (i.e., its Burgers vector is half of a 〈110〉 superdislocation), an anti-phase boundary (APB) forms behind it when the trailing partial enters γ' (Figure 7c). As the APB energy is much lower than the energy of the CSF, the APB widens, and the CSF remains narrow (Figure 7d), while the dislocation moves into the precipitate. In this process, the Co deposited on the γ-γ' interface ( Figure 7b,c) is drawn by the APB into the precipitate (Figure 7d).
The application example presented here demonstrates how the flexible structure of DAMASK allows to implement custom multiphysics models. In particular, it is shown how solute diffusion along dislocations ("pipe diffusion") and solute drag leads, under loading, to destabilization of the γ-γ' interface. In turn, this may give rise to precipitate dissolution, directional coarsening ("rafting"), loss of strength, and lifetime reduction. Therefore, the present efforts in alloy design are focused on the optimization of composition in a way as to delay or inhibit such effects.

Microstructure-Sensitive Prediction of Effective Properties
The macroscopic, average mechanical response of engineering alloys results from the interplay of the differently oriented grains and phases and their individual mechanical properties. Extensive testing is typically required to determine macroscopic properties such as orientation-dependent yield stress and ultimate tensile stress or plastic anistropy. Scale bridging or "virtual laboratory" approaches aim at partly replacing these experiments with microsctructure-sensitive simulation methods. For the determination of macroscopic properties, the exact representation of the microstructure is often computationally too expensive. Therefore, simplified representations are often used that allow to "only" determine effective properties without giving detailed insights into the deformation mechanisms. In the following, three simulation studies that predict macroscopic properties to exploit the full capacity of structural materials in engineering applications are presented.

A Two-Scale Scheme for Sheet-Metal Forming Simulations
The direct integration of microstructural details into componentscale simulations helps to improve their predictive quality. The simulation of sheet-metal forming operations, for example, requires to consider the plastic anisotropy and the spatially varying material behavior to predict phenomena known as "earing" and "springback", respectively. As direct simulations, where the microstructure is fully resolved, are prohibitive due to the extraordinary computational effort, a number of approaches have been proposed that introduce only selected microstructural aspects into sheet-metal forming simulations. "Texture-component crystal plasticity FEM" approaches [99,100] for example consider only the average crystallographic texture and neglect the spatial arrangement of the grains. Two-scale schemes can be either based on homogenization approaches, [101] on solving the mechanical boundary value problem on a "representative volume element" (RVE), or on model order reduction approaches [102,103] which can look up the response of the RVE from so-called offline-phase calculations.
Despite the high quality of microstructure-aware simulation approaches, using yield surfaces is still the de facto standard approach for sheet-metal forming simulations in industrial practice. A yield surface is defined in the 6D stress space as the boundary between purely elastic and elasto-plastic behavior. The von Mises yield criterion [104] and the Hill48 yield criterion [105] are widely used yield-surface descriptions for isotropic and anisotropic behavior, respectively. More involved yieldsurface descriptions improve the prediction quality but come at the cost of more elaborate testing protocols to adjust their higher number of parameters. Although this drawback can be tackled by (partly) replacing the experimental fitting procedure in "virtual laboratory" approaches with RVE simulations [106,107] similar to the offline-phase in model order reduction approaches, the description of the yield-surface evolution would still follow simple, phenomenological evolution laws developed by Hollomon, Swift, Voce, Ghosh, and Hockett-Sherby. Although these evolution laws allow to express the yield-surface expansion due to isotropic hardening and its translation in stress space due to kinematic hardening, shape changes associated with texture evolution require explicit consideration of lattice rotations.
It was the idea of Gawad et al. [108] to perform a yield-surface update in parallel to component scale simulations. In the  [60] www.advancedsciencenews.com www.aem-journal.com so-called hierachical multiscale scheme (HMS) scheme, the facet yield surface [109] is updated using the ALAMEL homogenization approach to predict the response of polycrystalline materials. An extension of this approach, where full-field crystal plasticity simulations replace the ALAMEL solution is presented in studies by Han et al. [110] During macroscale forming simulation, virtual experiments are conducted in an adaptive way to capture the plastic anisotropy evolution and update the yield function parameters accordingly.
Here, the HMS concept is demonstrated on the example of FEM simulations of a deep drawing test for an 2090-T3 aluminum alloy. As outlined in the studies by Han et al., [110] the mechanical response at the component scale is given by either the Yld2000-2d [111] or the Yld2004-18p [112] yield function. Before starting the simulation, the initial yield surface of the sheet alloy was calibrated based on virtual tests using the efficient spectral solver of DAMASK. [107,113,114] Then for each integration point of the macroscopic finite element model, the same DAMASK-based virtual tests are conducted to update the texture and corresponding yield surface at the material point individually at every 0.05 plastic strain increment. Thus the macroscopic FEM simulations take the evolution of the plastic anisotropy into consideration. In the Yld2000-2d yield function calibration, virtual tensile tests along 0 , 45 , and 90 relative to the rolling direction (RD) as well as equi-biaxial tensile tests are conducted to get eight data points (σ 0 , σ 45 , σ 90 , σ b , r 0 , r 45 , r 90 , and r b ). Eight equations can be constructed to obtain the eight parameters (α i , i ¼ 1…8) for Yld2000-2d simultaneously using a Newton-Raphson scheme. Similar virtual tests are conducted to calibrate the Yld2004-18p yield function; the details of the procedure are given in the studies by Han et al. [110] Figure 8 shows the initial texture as a (111) pole figure of the sheet alloy used for the cup drawing simulation. It is obtained from the volume element discretized by 10 Â 10 Â 10 ¼ 1000 grid points. The final configuration of the drawn cup simulated with the evolving Yld2000-2d yield function is shown in Figure 9. A strong deviation from the desired cup shape (i.e., "earing") Figure 9. Deformed configuration of the cup simulated with the evolving Yld2000-2d yield function started from the complete texture. The final textures at the three indicated positions are also shown as {111} pole figures, where the horizontal direction is aligned with the Transverse Direction (TD) and the vertical direction with the RD. Adopted with permission. [110] Copyright 2020, Elsevier. can be seen. The texture after deformation (represented by (111) pole figures) at three selected points is also presented in Figure 9. Significant and spatially heterogeneous texture evolution can be seen, indicating that the incorporation of texture is required to resolve the details of the cup profile. Figure 10 shows the evolution of the plastic anisotropy (represented by the uniaxial yield stresses, r-values, and the yield locus predicted by the Yld2000-2d yield function) at point P3 of the cup shown in Figure 9. It can be seen that the plastic anisotropy changes gradually as the plastic strain increases, and the change becomes significant at large strains. The scheme allowing the yield surface to evolve is a true multiscale modeling framework that can be applied to arbitrarily complex forming processes. It essentially combines the concepts of yield surface, virtual laboratory, and direct two scale approaches such as FEM 2 or FEM-FFT into an efficient and flexible microstructure-aware tool for simulating metal-forming operations. More details about the implementation can be found in the studies by Han et al. [110]

Assessment of Forming Limit Diagrams
For lightweight design of structural components, it is of utmost importance to exploit the full potential of a given material on the one hand but do not risk to deform the material in such a way that it damages. This requires a precise knowledge of the strains that a material can accommodate. A common approach to access the limits of a material in dependence of the applied principle strains (ε major and ε minor ) is the forming limit diagram (FLD) or forming limit curve (FLC). The critical strain, which separates the strain states that a material can accommodate from those that it cannot accommodate, is not the strain at fracture, but the strain at which localized necking occurs. The strain for fracture is usually higher than the strain for localized necking and is mainly governed by ε minor . [115] However, when the strain state is close to equi-biaxial tension, that is, ε major ¼ ε minor , fracture might occur without any localized necking. [116] Since significant efforts are associated with the experimental determination of FLDs, it is beneficial to replace them with predictive simulations.
Here it is shown-following [117] -how crystal plasticity simulations can be used to determine FLCs. To this end, DAMASK is coupled to the commercial FEM software LS-DYNA [118] to study the behavior of sheet metals under biaxial tension. Sample dimensions were selected as 4 mm Â 4 mm Â 1 mm and 8 mm Â 8 mm Â 1 mm (width Â breadth Â thickness). The geometry is discretized with linear brick elements with full integration capacity (eight integration points). An element size of 0.2 mm Â 0.2 mm Â 0.2 mm was selected such that the total number of elements is 20 Â 20 Â 5 ¼ 2000 and 40 Â 40 Â 5 ¼ 8000, respectively. Each integration point gets an individual orientation assigned and periodic boundary conditions are applied. Material parameters for aluminum listed in Table A.3, Supporting Information, are taken from the studies of Eisenlohr et al. [113] To determine the onset of necking, localized plastic deformation is quantified from a statistical point of view. The underlying procedure is shown in Figure 11 which shows the evolution of the cumulative probability functions of ε 11 for a random texture at plane strain tension condition, e.g., ε minor ¼ 0. Initially, the mean value, ε mean 11 , is approximately the same as the median value, ε median 11 . This indicates that ε 11 approximately follows a normal distribution. As the deformation progresses, ε mean 11 and ε median 11 deviate from each other and the distribution becomes more and more asymmetric. Eventually, ε mean 11 is mainly driven by the population above ε median 11 . These observations are consistent with the formation of a localized deformation in the direction of ε 11 (see the contour plots of ε 11 accompanying Figure 11). Based on the observations in Figure 11, the ratios ε mean ii =ε median ii (i ¼ 1 or 2), are used as the measure for the localization. As shown in Figure 11, ε mean ii =ε median ii ¼ 1 indicates (approximately) homogeneous deformation, and a deviation from 1 implies the onset of localization. However, from the contour plots in Figure 11, it can be seen that the transition from homogeneous deformation to localized deformation is a continuous process, which is consistent with the experimental observations made by Volk and Hora. [119] This makes it difficult to define the critical strain in an unambiguous way. Figure 12 shows ε mean 11 =ε median 11 and the macroscopic stress, P ij , computed as the average over all the integration points in dependence of the applied strain. This allows to compare the statistical Figure 10. Evolution of the plastic anisotropy at point P3 of the cup simulated with the evolving Yld2000-2d yield function started from the complete texture. a) Normalized yield stress versus angle to RD Θ. b) r-value versus angle to RD Θ. c) Yield locus. The symbols are the data points obtained from the corresponding virtual experiments. Adopted with permission. [110] Copyright 2020, Elsevier.
www.advancedsciencenews.com www.aem-journal.com criterion to Drucker 's stability criterion, [120][121][122][123] which states that the deformation instability occurs when strain softening takes place, that is, the forming limit is at the maximum stress on the stress-strain curve. As observed in Figure 12, the maximum stresses do not occur exactly when ε mean ii starts to deviate from ε median ii , but when ε mean ii =ε median ii ¼ 1.05 $ 1.10. This implies that the strain hardening capacity in the localized volume is not fully exhausted immediately at the initiation of the localization, but can be sustained up to a certain extent.
The simulated mean and extreme values for FLDs of six typical fcc texture components (see studies by Ma et al. [124] for details on the naming scheme) are shown in Figure 13. Five statistically equivalent orientation distributions are used for each texture component and the FLDs are determined based on Drucker's criterion. A strong and systematic influence of the texture on the predicted FLDs can be seen. However, the calculated FLDs still exhibit a noticeable scatter. The scatter introduced by the different orientation distributions is not surprising as the localized deformation depends highly on the local spatial orientation distribution. Figure 13 also shows that the FLDs from the 40 Â 40 Â 5 element mesh are slightly lower than those with 20 Â 20 Â 5 elements. A possible reason is the use of periodic boundary conditions. If a localized deformation band is formed, the orientation or the shape of the band should follow the imposed periodicity. Thus, periodic boundary conditions impose restrictions to the localized deformation. In an RVE with more degrees of freedom, there is more freedom for the localized deformation band to adopt its orientation or shape as it is less restricted by the imposed boundary conditions.
The texture dependence revealed in Figure 13 is that 1) the cube texture component has the best formability, and it is the only orientation distribution whose forming limit increases when ε major increases; 2) random, copper, and S have comparable formability to each other, and all of them are better than Goss and Brass; 3) Goss has the poorest formability.  versus ε major (top) and the macroscopic stresses P ij versus ε major (bottom) at plane strain tension load (ρ ¼ 0.0). The solid circles mark the maximum stresses. Each line corresponds to a random orientation distribution instance and all simulations were performed using the 20 Â 20 Â 5 elements mesh.
www.advancedsciencenews.com www.aem-journal.com The texture dependence as shown in Figure 13 is actually in good agreement with Figure 4f in ref. [125]. This agreement is surprising as the origin of the localized deformation in the study presented here is the microstructural heterogeneity, whereas in ref. [125], it is a predefined geometrical imperfection. It can be therefore concluded that the governing effect is only the crystallographic texture. Therefore, crystal plasticity simulation can be considered as an suitable alternative to experimental tests for the determination of FLDs.

Microstructure Influence on r-Values and Hardening of Dual-Phase Steels
Dual phase (DP) steels, consisting of a soft ferritic and a hard martensitic phase, are an important family of steel grades used by the automotive industry. [126] The ferrite is mainly responsible for the ductility, whereas the martensite give these steels their strength. Typical DP microstructures comprise of complexshaped ferrite grains surrounded by smaller martensite particles. The mechanical properties of the microstructure depend on the volume fraction of the phases, the size distributions of grains and particles, their spatial distributions, as well as the crystal orientations of the constituents. One of the important mechanical characteristics for deep-drawing applications is the r-value or Lankford coefficient. [127] An accurate prediction of hardening behavior and r-value is required for faster product development of improved DP steel grades.
To create RVE that incorporate the complex features of DP steels, Tata Steel has developed a microstructure generation tool which is based on the multi-level Voronoi tessellation [128] and has been enriched with additional features. In contrast to the standard Voronoi tessellation, which is only capable of creating convex shaped cells (i.e., grains), the multi-level algorithm makes it possible to create complex-shaped grains. Moreover, different size and spatial distributions for ferrite grains and martensite inclusions can be generated. More details of the tool can be found in refs. [129,130].
In the study presented here, the influence of the spatial distribution of martensite particles and their orientations on hardening behavior and the r-value is investigated. 2D electron backscatter diffraction (EBSD) data from an industrial DP600 steel grade forms the basis for building 3D microstructures. Figure 14 shows the measured microstructure colored according to the inverse pole and the ferrite grain size distribution. The area fraction of martensite was determined as 11.7%. The orientation distributions of ferrite and martensite are shown in Figure 15. For the RVE generation, the phase volume fraction, ferrite and martensite orientation distribution, ferrite and martensite size distribution, as well as the spatial distribution of martensite has been used. To determine the volume fraction of martensite from the area fraction the simple stereological formula is used. In this, dV i and dA i are respectively the volume fraction and area fraction of phase i, and V and A are respectively the total www.advancedsciencenews.com www.aem-journal.com volume of the RVE and the total area of the EBSD scan. A realistic texture can be sampled from 2D EBSD data on basis of the wellknown texture theory of Bunge. [131] To this end, standard cleaning and thresholding operations available in the EBSD software are used to group the original pixels into grains. Every grain as well as every martensite particle has been assigned a single orientation. The grain size and the discrete orientation data were then used to assign orientations with the proper intensities to the grains in the 3D RVE. According to ref. [131], a texture can be described with an orientation distribution function (ODF) f(g), which is defined as where dV is the fraction of sample volume with orientation g in orientation bin dg. V is the total sample volume and φ 1 , Φ, φ 2 are the Euler angles. Reinterpreting dA i as the area of grain i and dV i as the volume of all grains having the same orientation as grain i allows to use Equation (1) also for determining the volume distribution in the 3D microstructure. A Monte Carlo optimization procedure was used to map orientations to grains such that the volume fraction of each orientation comes as close as possible to the area fraction of the same orientation in the EBSD scan. Figure 16 shows three different RVEs for a DP steel grade with exactly the same volume fraction of martensite: 11.7%. All RVEs have been built with 72 000 Voronoi cells which are grouped into 302 complex-shaped grains in a cube of 50 Â 50 Â 50 μm. Moreover, in cross sections, an apparent size ratio of 6:1 between ferrite and martensite is observed, which is in agreement with the EBSD measurement (the experimentally measured average ferrite grain diameter is 5.9 μm and the average martensite diameter is 1.0 μm). This requires to have an average grain diameter of 8.5 μm for ferrite and of 1.5 μm for martensite in the 3D RVE.
To study the influence of the crystallographic orientation of martensite and ferrite, three different distributions have been generated: a) Random orientations of ferrite grains and martensite particles; b) Ferrite having the same texture as measured by EBSD and martensite particles being assigned the orientation of a neighboring ferrite grain; c) Ferrite grains as well as martensite particles having identical texture as measured by EBSD.
To study the influence of the spatial distribution of the martensite particles, three different cases have been studied: 1) Martensite particles are distributed partly in bands and partly random along ferrite grain boundaries. This resembles best the experimental microstructure in the EBSD measurement; 2) A fully banded distribution of martensite; 3) A completely random distribution of martensite. The microstructures are discretized into a regular grid for use with the FFT-based spectral solver of DAMASK. [19,113,114] For both phases the phenomenological crystal plasticity model is used. Constitutive parameters for ferrite-given in Table A.5, Supporting Information-are obtained from a nanoindentationbased optimization procedure for a forming simulation. The most relevant crystal plasticity parameters for martensite, namely the initial hardening h 0 and resistances ξ 0 and ξ ∞ have been adjusted such that holds. Here, σ sim ðε i Þ and σ exp ðε i Þ are the simulated and the experimentally obtained stress response for a given strain ε i under uniaxial loading. The resulting parameters are given in Table A.6, Supporting Information. The mechanical properties of martensite are, however, strongly influenced by its carbon content. A higher carbon content will make the martensite harder. Rodriguez and Gutierrez [132] explain how hardening curves of both phases can be related to their chemical composition. To estimate the carbon content of martensite in dependence of its volume fraction, data from the experiments of Thomser et al. [133] have been used here. This allows to derive a relation between volume fraction of martensite, its carbon content, and the crystal plasticity parameters h 0 , ξ 0 , and ξ inf . In the considered range, a linear relation for all three parameters has been found. This procedure allows to perform simulations over a broad range of compositions without the necessity to perform additional tests for parameter identification. Figure 17 shows the hardening behavior in dependence of the assigned texture and the microstructure model. A good agreement with the experimental data is achieved for the most realistic case (model 1) independently of the texture. Although a complete banding (model 2) gives a significant higher stress response, the response of a spatially random distribution (model 3) gives only a slightly softer response as the realistic case (see Figure 17c). As expected, the loading angle (0 and 45 with respect to the RD) has no influence if a random texture is used (see Figure 17a). This changes if a realistic texture is assigned to the RVE. Comparing Figure 17b with 17c clearly shows that especially the martensite orientation has a significant influence on the hardening rate. With a realistic texture, the hardening response for loading perpendicular to the RD is higher than for loading along the RD for all the models. Figure 18 shows the r-value for loading angles of 0 , 45 , and 90 with respect to the RD in dependence of the assigned texture and the microstructure model. Experimentally obtained values are also provided. It can be seen that the spatial distribution of the martensite has a less pronounced effect than the crystallographic orientation. Although no systematic trend between the random orientation distribution (Figure 18a) and the realistic texture for ferrite (Figure 18b) can be seen, using also the measured orientations for martensite (Figure 18c) results in a relevant increase in the r-values. While the simulated r-values are too high when using the measured texture, the trend for the realistic microstructure (model 1) is in almost perfect agreement with the measurement.
The use case studied here shows how advanced microstructure generation tools that enable the incorporation of experimentally obtained microstructure descriptors in conjunction with crystal plasticity simulations can substitute the mechanical tests. This is of special importance for materials with complex microstructures such as DP steels as it allows to investigate systematically the influence of individual microstructural features on various mechanical properties. www.advancedsciencenews.com www.aem-journal.com

Conclusion and Outlook
The development of new alloys with improved properties and the economically efficient use of structural materials with small safety margins require the precise knowledge of the materials behavior in dependence of their microstructure. As the individual microstructural constituents such as second-phase particles, grains, subgrains, and dislocations respond differently to the thermo-chemo-mechanical loading to which the material is subjected during production and in service, combined modeling of these diverse effects is required for ICME. In this overview article, it was demonstrated how multiphysicsenhanced crystal plasticity simulations can be used to study the mechanical behavior of metallic materials at different length scales. All simulations have been performed with DAMASK, a free and open-source simulation toolkit. Investigations at the grain or subgrain scale benefit especially from the inclusion of additional physical effects. In Section 2.1 and 2.2, it is presented how a coupled phase field scheme for fracture and twinning, respectively, enables a very detailed perspective on the microstructure evolution under mechanical loads. This allows the development of microstructures with second-phase particle distributions which have an increased strength and damage resistance. As outlined in Section 2.3, the structure of DAMASK can be also used to implement simulation approaches beyond the original scope of grain scale crystal plasticity. More specifically, an APFCM has been implemented and was used to study the interaction of solutes with dislocations and stacking faults. For the determination of effective properties, such multiphysics enhancements are often not required and are disabled to increase the computational performance. As shown in Section 3.2 and 3.3, plain crystal plasticity simulations are sufficiently accurate to predict the macroscopic properties such as the r-value and the forming limit based on crystallographic texture and microstructure. Hence, they are a cost-effective alternative to experimental investigations when evaluating the performance of existing materials. The obtained properties can be directly used in engineering FEM to design and evaluate components made out of the probed material. As shown in Section 3.1, an even closer coupling between component-scale simulations and material mechanics is feasible and beneficial for improving the prediction quality of the large scale simulations.
The aforementioned examples clearly demonstrate how ICME tools in general and DAMASK in particular can be used to study diverse materials science and engineering related questions. Especially noteworthy is the success in combining different models into holistic, general purpose simulation tools. Nevertheless, having the complex microstructures in modern high-performance materials and the various thermo-chemo-mechanical loads in mind, it becomes clear that there is still space for significant improvements. Integration of models for recrystallization and grain growth, for example, is an essential prerequisite for through-process modeling.