Loading next page...
/lp/arxiv-cornell-university/computational-homogenization-of-polycrystals-pBLJ97tIIW
- ISSN
- 0065-2156
- eISSN
- ARCH-3331
- DOI
- 10.1016/bs.aams.2018.07.001
- Publisher site
- See Article on Publisher Site
Abstract
This paper reviews the current state-of-the-art in the simulation of the mechanical behavior of polycrystalline materials by means of computational homogenization. The key ingredients of this modelling strategy are presented in detail starting with the parameters needed to describe polycrystalline mi- crostructures and the digital representation of such microstructures in a suit- able format to perform computational homogenization. The di erent crystal plasticity frameworks that can describe the physical mechanisms of deforma- tion in single crystals (dislocation slip and twinning) at the microscopic level are presented next. This is followed by the description of computational homogenization methods based on mean- eld approximations by means of the viscoplastic self-consistent approach, or on the full- eld simulation of the mechanical response of a representative polycrystalline volume element by means of the nite element method or the fast Fourier transform-based method. Multiscale frameworks based on the combination of mean- eld ho- mogenization and the nite element method are presented next to model the plastic deformation of polycrystalline specimens of arbitrary geometry under complex mechanical loading. Examples of application to predict the strength, fatigue life, damage, and texture evolution under di erent conditions are pre- sented to illustrate the capabilities of the di erent models. Finally, current challenges and future research directions in this eld are summarized. Keywords: Homogenization theory; crystalline solids; multiscale modelling; crystal plasticity. Corresponding author Preprint submitted to Advances in Applied Mechanics July 24, 2018 arXiv:1804.02538v2 [cond-mat.mtrl-sci] 22 Jul 2018 2 1. Introduction Physically-based modelling of the mechanical behavior of materials is an integral part of Materials Science and Engineering, as an indispensable tool to understand the relationship between microstructure and properties. However, the complexity of the deformation mechanisms that need to be represented, in addition to limitations in computational resources, have un- til recently restricted the goal of these formulations to provide qualitative trends or to complement phenomenological models that need to be adjusted from experimental data obtained under di erent conditions. Recent advances in simulation techniques, computational power, and multiscale modelling strategies are rapidly overcoming these limitations. These new paradigms are aimed at developing formulations that can quantitatively predict rela- tionships between processing, microstructure and properties, and contribute to the design of materials with a given set of properties in silico, before they are manufactured in the laboratory. Successful examples of the latter can be found in the realm of metallic alloys (Allison et al., 2006; Olson, 2013) or structural composites (LLorca et al., 2011; Gonz alez et al., 2017), and the roadmaps for implementation of these strategies in di erent industrial sectors are clearly delineated (Allison et al., 2013). Although the number of success stories is still limited, they show a dramatic reduction in the time necessary to develop new materials with optimized properties for speci c applications, as well as a large return-of-investment. Thus, there is a manifest interest in both the scienti c community and industry to improve and develop new multiscale modelling strategies. Within this scenario, homogenization theory stands out as one of the most important tools to relate microstructure with e ective properties in heterogeneous materials (Nemat-Nasser and Hori, 1993; Torquato, 2001). Homogenization theory assumes that the characteristic length scale of the microscopic domains (the average size of the heterogeneities) is much larger than the molecular dimensions (and, thus, continuum mechanics is applica- ble), and much smaller than the characteristic length scale of the macroscopic sample. Under these conditions, the macroscopic or e ective properties of the material can be determined from geometrical features of the microstruc- ture and properties of the di erent heterogeneities in the material. Following the pioneer work of Eshelby (1957), di erent homogenization mean- eld for- mulations were developed rst for linear materials (Kr oner, 1958; Hashin and Shtrikman, 1963; Hill, 1965b; Mori and Tanaka, 1973). These methods 3 were later extended to deal with non-linear mechanical behaviors (nonlinear elasticity, plasticity and viscoplasticity) by means of di erent linearization approximations (Hill, 1965a; Berveiller and Zaoui, 1979; Tandon and Weng, 1988; Ponte-Castaneda, ~ 1991, 1996). Among these mean- eld approaches, the elastic self-consistent (SC) scheme (Hershley, 1954) arises as the most reliable approach to obtain estimations of the elastic response of polycrystals (Hill, 1965b). This approximation assumes that each domain in a set that statistically represents the heteroge- neous material is embedded in a homogenous e ective medium whose elastic properties are not known a priori, but need to be obtained as an average over those heterogeneities. The linear elastic SC formulation was extended to non- linear viscoplastic polycrystals by Molinari et al. (1987a), taking into account that plastic deformation in each grain is accommodated by dislocation glide, resulting in local mechanical responses determined by the crystallography and lattice orientation of each grain. This methodology provides a robust framework to establish physically-based relationships between texture, mi- crostructure and mechanical properties of plastically deforming polycrystals, accounting for their evolving microstructure and the anisotropy of the single crystal grains. The viscoplastic self-consistent (VPSC) approximation was later improved by Lebensohn and Tom e (1993); Lebensohn et al. (2007), and the VPSC formulation is nowadays widely used to simulate the mechanical behavior of polycrystalline aggregates. However, homogenization strategies based on mean- eld approximations have two main limitations. The rst one is that these approaches are based on a description of the microstructure based on average values of grain sizes, shapes and orientations, and therefore cannot take into account the e ects of local heterogeneities such as clustering of grains with di erent size, shape, or orientation, or to represent a given misorientation distribution. While these local microstructure features has often negligible e ects on properties that depend on the average values of the micromechanical elds (like, for instance, the elastic modulus), this is not the case for properties that de- pend on the extremal values (e.g. damage) (Segurado et al., 2003; LLorca and Segurado, 2004; Segurado and LLorca, 2005). The second limitation arises from assuming that the micromechanical elds are constant in each phase/grain. Under these circumstances, it is dicult to simulate phenom- ena in which these elds are localized in narrow bands in one or several phases/grains, as a result of plastic deformation, damage progression, low strain rate sensitivity (i.e. high non-linearity), strong mechanical contrasts, 4 etc. (Gonz alez et al., 2004; Gonz alez and LLorca, 2007; Totry et al., 2008). More sophisticated non-linear homogenization methods have been developed in recent years, using linearization schemes at phase/grain level that also in- corporate information on the second moments of the eld uctuations in each phase/grain, to ameliorate this second limitation, at the cost of more com- plex and expensive numerical algorithms (de Botton and Ponte-Castaneda, ~ 1995; Liu and Ponte-Castaneda, ~ 2004). These limitations of mean- eld homogenization approaches can be over- come by means of full- eld homogenization. Using this strategy, the e ec- tive properties of the polycrystal are determined by means of the full- eld solution of a boundary value problem of a microstructural Representative Volume Element (RVE) under homogeneous boundary conditions. It has been established that the success and accuracy of computational homoge- nization approaches for polycrystals rely on the adequacy of three factors: the representation of the microstructure, the constitutive description of the single crystal behavior, and the numerical approach to solve the boundary value problem. From the RVE perspective, the rst requirement is that the size of the RVE has to be large enough to provide an accurate statistical representation of the polycrystal's microstructure, as well as to lead to e ec- tive properties that are independent of the RVE dimensions. Obviously, the larger the RVE and the more complex the single crystal constitutive behav- ior, the higher the computational cost and, thus, fast and ecient numerical procedures to solve the boundary value problem are essential. Full- eld homogenization has been traditionally carried out using the - nite element method (FEM) to solve the governing partial-di erential equa- tions (PDEs) of micromechanics, in combination with the CP formalism de- veloped by Peirce et al. (1982a). This model takes into account the process of plastic deformation by crystallographic slip and provides an accurate and physically-based representation of this phenomenon at the microscopic level within each grain (Becker, 1991; Miehe et al., 1999; Delaire et al., 2000; Raabe et al., 2001; Bhattacharyya et al., 2001). Alternatively, Moulinec and Suquet (1994) proposed a spectral formulation based on the ecient Fast Fourier Transform (FFT) algorithm to perform full- eld homogenization of a periodic RVE. The governing PDEs were solved in this case by means of the computation of convolution integrals involving the periodic Green's func- tion associated with the displacement (or velocity) eld of an homogeneous linear reference material. This strategy was initially applied to composite materials (Moulinec and Suquet, 1998; Eyer and Milton, 1999; Michel et al., 5 2000), in which the source of heterogeneity is related to the spatial distribu- tion of phases with di erent mechanical properties. It was later adapted to polycrystals, in which the heterogeneity is related to the spatial distribution of crystals with anisotropic mechanical properties (Lebensohn, 2001; Leben- sohn et al., 2005, 2008, 2009). It should also be noted that grain boundaries (GBs) of very di erent character are also ubiquitous in polycrystals and have to be modelled appropriately. As a consequence of the above developments, simulation of the mechani- cal behavior of polycrystals has signi cantly improved in the last 3 decades and mature strategies based on either mean- eld or full- eld homogenization approaches are nowadays available. The state-of-the-art of Computational Homogenization of Polycrystals (CHP) is summarized in this review, which is structured as follows. Section 2 is focused on the description of polycrys- talline microstructures, including the experimental techniques available for their determination, and the digital representation of microstructures in a suitable format to perform the computational homogenization. Section 3 presents the Crystal Plasticity (CP) framework used by the di erent homog- enization methods to describe the plastic deformation of each single crystal within the polycrystal. Both phenomenological and physically-based models are described, as both are of interest for di erent applications. A detailed description of the Green's function-based VPSC homogenization approach is given in Section 4, while full- eld homogenization strategies based on the numerical solution of a boundary value problem of an RVE using either the FEM or FFTs are presented in sections 5.1 and 5.2, respectively. While a detailed description of CPFEM can be found elsewhere (Roters et al., 2010b,a) and, therefore, is omitted in this review, the ecient FFT-based formulation is described in detail. In particular, the underlying commonal- ities and di erences between both Green's function-based CP formulations (VPSC and FFT-based models) are presented with an uni ed notation. Fur- thermore, section 6 presents the formulation of multiscale frameworks for modeling the plastic deformation of polycrystalline specimens of arbitrary geometry under complex mechanical loading. These frameworks are based on the combination of polycrystal homogenization and FEM, by connecting the single crystal (microscopic) behavior with the specimen (macroscopic) response through an intermediate (mesoscopic) scale. The behavior at the mesoscopic scale is obtained from the homogenized response of a polycrystal, that represents each (polycrystalline) material point in the FE mesh. A few applications of the described models are presented in Section 7, showing the 6 capabilities of the di erent simulation strategies. These examples illustrate virtual design strategies, in which the mechanical response of polycrystalline microstructures can be accurately computed and optimized in silico before the material is actually manufactured. Finally, potential new applications of CHP are summarized in the last section, together with current limitations and topics that should be addressed in the future. Throughout this paper, vectors are denoted by bold lowercase roman letters a, second-rank tensors by bold capital roman letters or bold Greek letters (A or ) and fourth-rank tensors by A. In addition, bold capital roman (E) and greek letters ( ) are used to express volume-average vectors, such as stress, strain or velocities. A Cartesian coordinate system is used with respect to the orthonormal basis (e ; e ; e ). The notations for tensor 1 2 3 product, contraction and double contraction products are: a b = a b e i j i e ; AB = A B (e e ); a b = a b and A : B =A B . Finally 1 and I j ik kj i j i i ij ij stand for the second- and fourth-rank identity tensors, respectively. 2. Representation of the microstructure The rst step to perform computational homogenization is the genera- tion of synthetic RVEs of the microstructure. The de nition of an RVE can be approached from two di erent perspectives. On the one hand, it can be de ned as the minimum volume of the microstructure whose properties (obtained using either homogeneous traction or displacement boundary con- ditions) are equal to the e ective properties of the heterogeneous solid. On the other hand, the RVE can be de ned as a the minimum volume of the microstructure whose statistical descriptors (that de ne the features of the heterogeneous microstructure) are equivalent to those found in the heteroge- neous solid. Another important concept is that of Statistical Representative Volume Element (SRVE) (Jeulin et al., 2004). The SRVE does not ful ll one or both conditions imposed to the RVE by itself, but it does ful ll them in a statistical sense: the e ective properties and statistical descriptors of the mi- crostructure obtained by averaging the results obtained for many SRVE are equivalent to those obtained for an RVE. The use of SRVE is sometimes very useful in computational homogenization because the simulation of very large RVEs is computationally much more expensive than that of many smaller SRVEs. 7 2.1. Microstructural parameters Polycrystals are aggregates of single crystal grains, which are the building blocks of the RVE. Crystals are modelled as homogeneous polyhedra with well-de ned crystallographic orientations separated by grain boundaries. In addition, internal discontinuities in the lattice of a grain due to the presence of twinned regions can also be included, due to the importance of these fea- tures to reproduce the mechanical response at the local level (for instance, in the nucleation of fatigue cracks). The main features to generate an RVE of a polycrystalline microstructure are the grain size and shape distribution, together with the Orientation Distribution Function (ODF) that de nes the crystallographic texture of the polycrystal. In addition, RVEs can include local microstructure heterogeneities or speci c grain misorientation distribu- tions. Optical microscopy of polished sections of the microstructure is the stan- dard technique to measure the grain size and shape distribution. Chemical etching can be used to reveal the grain boundaries and the statistical param- eters that characterize the grain distribution can be automatically obtained using image-analysis software tools. Of course, the 2-D information obtained has to be transformed into 3-D and this operation may or may not be simple, depending on the microstructural features. For instance, there are algorithms available that can infer the 3-D grain size distribution from the data obtained in 2-D for equiaxed grain distributions (Heilbronner and Bruhn, 1998). Sim- ilarly, grain shapes can be accurately obtained from 2-D images in di erent orientations in the case of materials with well-de ned anisotropy, as it is the case of rolled and extruded sheets. From the viewpoint of grain orienta- tion, the ODF can be easily obtained from X-ray di raction using texture goniometers. The information obtained can be directly used to model grain orientation distribution in the RVE. More sophisticated methods have been developed in recent years to pro- vide a more accurate quanti cation of the microstructure in polycrystals. Au- tomated serial sectioning by means of a microtome (Alkemper and Voorhees, 2001) or mechanical milling (Spowart et al., 2003; Spowart, 2006) can be used in combination with software tools to build the 3-D microstructure and extract microstructural features, to determine grain size and shape distribu- tions in 3-D. If higher resolution is required, sequential milling can be carried out using a focus ion beam or a femtosecond laser within a scanning electron microscope (Uchic et al., 2006; Echlin et al., 2011). One additional advantage of this latter approach is that 3-D chemical and crystallographic information 8 can be obtained during the process. In the particular case of polycrystals, the use of Electron Back-Scattered Di raction (EBSD) provides information about the crystallographic orientation of each grain, the orientation relation- ship at the grain boundaries as well as the presence of twins within the grains (Uchic et al., 2004; Konrad et al., 2006; Fern andez et al., 2013). While serial sectioning techniques are destructive, non-destructive characterization of the 3-D structure of polycrystals can be obtained by means X-ray tomography. Standard phase-contrast tomography provides information about porosity and second phases. Moreover, information about grain size and crystallo- graphic orientation is presently available by means of di raction-contrast tomography (Schmidt et al., 2008; Proudhon et al., 2015). Thus, detailed quanti cation of the microstructure of polycrystals is nowadays possible, but the acquisition time and the management and interpretation of the infor- mation contained in the massive datasets generated by 3-D characterization techniques may be problematic. Thus, the microstructural information nec- essary to build the RVE for computational homogenization should be assessed beforehand, to ensure ecient data collection and reconstruction. It should be nally noted that the relevant microstructural information can be captured from 2-D or 3-D datasets by statistical descriptors, such as the 1-point and 2-point statistics (Torquato, 2001). 1-point statistics cap- ture the probability associated with nding a speci c phase or orientation at a point thrown randomly into the microstructure. They can be used to characterize the volume fraction of di erent phases or the texture of the polycrystal. 2-point statistics indicate the probability that both ends of a segment of given length thrown randomly into the microstructure lie in the same phase. They can be used to characterize the spatial distribution, size and shape of the di erent grains in the polycrystal. Once the particular sta- tistical descriptors of a microstructure have been determined, phase-recovery algorithms can be used to generate RVEs, which are statistically equivalent, i.e. they share the same statistical descriptors (Fullwood et al., 2008; Niez- goda et al., 2010; Do sk a r et al., 2014). The use of 2-point correlation function is particularly interesting for describing composites or two phase-materials. In the case of polycrystals Niezgoda et al. (2010) each orientation should be considered as a di erent phase and the number of 2-point correlation func- tions to describe the microstructure becomes large. In this case, alternative microstructure descriptors such as the probability distribution of grain sizes, or grain principal axis distribution are commonly used. Motivated by this observation, we further explore the potential of Wang 9 tiles to represent long-range spatial correlations in disordered microstruc- tures, a problem common to materials science [11], geostatistics [12], or im- age analysis [13]. In this regard, two closely related applications can be dis- tinguished, namely, the microstructure reconstruction [14?16] based on given spatial statistics and microstructure compres- sion [17?19] aiming at ecient representation of material structures in multiscale computations 2.2. Digital representation of the microstructure Once all the critical information about the microstructure of the poly- crystal has been gathered, it is necessary to build a digital model of the RVE to perform computational homogenization. The RVEs can be constructed as a one-to-one representation of an actual microstructure measured from X-ray computed tomography or serial-sectioning data, or by generating synthetic microstructures from the statistical descriptors representing the microstruc- ture. Two di erent types of discretizations are usually employed to de ne the geometry of the RVE. In the rst one, the RVE is divided in voxels and grains are formed by groups of contiguous voxels that have the same crys- tallographic orientation (Fig. 1a). Voxel-based discretizations are the best option for direct representation of microstructures measured by three dimen- sional techniques because they can be directly extracted from the measured data. In addition, voxel-based discretizations can be directly exported into FFT-based codes as grid points or into FEM as regular hexaedral elements. Such voxel-based nite element meshes have good quality metrics and can be deformed up to very large strains (Segurado and LLorca, 2013). These voxel-based discretizations present, however, two major drawbacks. Firstly, the grain boundaries are stepped surfaces and, thus, this type of representa- tion is not appropriate to simulate phenomena localized at GBs, like grain boundary sliding. Secondly, the voxel-based FEM discretization often leads to very large number of elements because the voxel size is controlled by the dimensions of the smallest features that have to be resolved within the RVE. These limitations of the voxel-based digital models can be overcome with representations of the grain structure based on tessellation. A tessellation is a subdivision of the 3-D space into convex polyhedra that intersect only at their boundaries, which are at surfaces. The most popular one is the Voronoi tessellation, which is generated by a set S of points in the three dimensional space by assigning a volume V to the point P (x ) 2 S formed x i i by all points P (y) in the space which have P as their nearest neighbour (Redenbach et al., 2012). Cell boundaries in the Voronoi construction are, 10 Figure 1: Digital representations of an RVE of a polycrystal. (a) Voxel-based representa- tion. (b) Voronoi-based representation. however, always equidistant from the generators of their cells and the range of cell patterns which can be generated is limited. Thus, Voronoi tessella- tions are not always suitable to reproduce the actual grain size distribution and weighted generalizations of the Voronoi model are frequently used. One possible generalization is the Laguerre tessellation in which the volume V in the space associated to the point P 2 S is formed by the points P (y) that full ll the condition P (y) 2 V if d (y; x ) < d (y; x ); j 6= i and x 2 S (1) x L i L j j where d (y; x ) is the "Laguerre" distance between points y and x , which L i i is given by 2 2 d (y; x ) = kx yk r (2) L i i where r (> 0) is the weight associated to point P . This de nition leads to a i i partition of the space formed by convex, space- lling polyhedrons. It can be shown (Lautensack, 2007; Xue et al., 1997) that if S is chosen as a system of nonoverlapping spheres characterized by the coordinates of their centers, x and the corresponding radius, r , each cell of Laguerre tesselation completely 11 contains its generating sphere and the volume distribution of Laguerre cells is almost equal to volume distribution of their generating spheres. Thus, the strategy to generate grain structures within the RVE begins with the experimental grain size distribution, which can be often approxi- mated by a log-normal function. A polydisperse sphere distribution, follow- ing the experimental grain size distribution, is then introduced in the RVE and densely packed using collective rearrangement algorithms (Torquato, 2001) such as the force biased algorithm (Bargie l and Mo scinski, 1991). This algorithm starts with an initial distribution of spheres S(x ; r ) character- i i ized by the position of the center x and the radius r distributed in the i i RVE. In this stage, overlapping of spheres is possible and allowed. Then, the algorithm attempts to reduce the overlaps between spheres by pushing apart overlapped spheres while small spheres are pushed to ll the empty spaces between large ones. After certain number of iterations, repositioning of overlapped spheres is stopped and the spheres gradually shrank to reduce the total amount of overlaps below a certain threshold. Finally, the coordi- nates of the centers of the spheres and their diameter are provided as output for the Laguerre tesselation. In other cases, Monte Carlo algorithms were used to obtain the spatial distribution of the set points for the tessellation, so the nal cell size distribution coincides with the experimental grain size distribution (Cruzado et al., 2015; Mandal et al., 2018). Once the grain size distribution has been reproduced in the RVE, the crystallographic texture can be introduced by assigning di erent orientations to the grains to reproduce the statistical distribution given by the ODF. It should be noted that the minimum number of grains in the RVE should be large enough to reproduce accurately both the grain size distribution and the texture of the polycrystal. Moreover, more sophisticated grain orientation algorithms can be used to account for the presence of a given fraction of low- or high-angle grain boundaries, which may lead to important di erences in the mechanical behavior at both microscopic and macroscopic levels. The digital representation of the microstructure is an important and time consuming task and microstructure builders have been developed, such as Neper (Neper, 2018) and Dream3D (Dream.3D, 2018). They also provide tools to clean up the voxelized microstructures obtained from tomography or serial-sectioning or from tessellation and to discretize the microstructure for full- eld simulations. 12 3. Crystal plasticity models The rst model that described the plastic deformation of metallic single crystals as a result of crystallographic slip was proposed by Taylor and Elam (1923, 1925). In this seminal work, the deformation of Al single crystals was analysed and explained as the result of the shear deformation along twelve slip systems and the driving force for the shear deformation was the resolved shear stress on each slip system. A few years later, this model was used as to analyze the deformation of a polycrystal as an aggregate of grains (Taylor, 1938). The ideas of Taylor for the deformation of single crystals were adapted into the framework of continuum mechanics by Hill (1966), in the case of small strains. The theory, based on the general internal variable thermodynamic formalism, was extended to nite deformations in the 70s by Rice (1971) and Hill and Rice (1972). They used the concept of the multiplicative decomposition of the deformation gradient into elastic and plastic parts, introduced by Lee and Liu (1967). Di erent constitutive models based on this framework were developed for single crystals in the 80's using either rate-independent formulations (Peirce et al., 1982b; Peirce, 1983) or viscoplasticity (Asaro and Needleman, 1985). In parallel, attention was also paid to the development of rigorous numerical implementations of the models, including ecient and well-posed integration methods for the highly nonlinear viscoplastic laws (Cuitino ~ and Ortiz, 1992) and rigorous integration for the nite deformation framework (Miehe, 1996). The result of all these studies { and many more not reviewed here { is a well-stablished theory of CP which will be summarized below. The starting point of most CP models is the multiplicative decomposition e p of the deformation gradient F into its elastic (F ) and plastic (F ) parts (Fig. 3) e p F = F F (3) where the con guration de ned by F is called the relaxed or intermediate con guration. In the context of CP, it is assumed that F leaves the crystal lattice undistorted and unrotated (Rice, 1971; Hill and Rice, 1972) and the rotation of the lattice is determined only by F . Although this decomposition is generally accepted for CP, several issues as the existence and uniqueness of the decomposition and its connection with the microscopic distortion gener- ated by the dislocations are still under debate (Reina and Conti, 2014; Reina et al., 2016). 13 𝛼 𝑭 𝑭 𝒎 * 𝒆 𝒑 𝑭 = 𝑭 𝑭 Figure 2: Schematic of the multiplicative decomposition of the deformation gradient From the de nition of the velocity gradient, L, expression (3) leads to 1 1 1 1 e e e p p e _ _ _ L = FF = F F + F F F F (4) p p p where L = F F stands for the plastic deformation rate in the intermedi- ate con guration. The constitutive equations can be obtained from the energy density per unit volume expressed in the intermediate con guration (Cuitino ~ and Ortiz, 1992; Han et al., 2005a). Following the internal variable formalism, the free energy density, , can be written as = (F; F ; q) (5) where q is a set of internal variables. The free energy density can be split into the elastic and plastic energy densities (Cuitino ~ and Ortiz, 1992), e p p p = (FF ) + (F ; q): (6) and the second Piola-Kircho stress tensor in the intermediate con guration S can then be obtained as S = (7) @E where E stands for the Green-Lagrange elastic strain in the intermediate con guration and is given by 1 T e e e E = (F F I): (8) 14 The second Piola-Kircho stress in the intermediate con guration S is related with the Cauchy stress, , according to 1 1 T e e S = F F (9) where J is the determinant of F. Assuming a quadratic potential for the elastic energy in eq. (6), S can be expressed as a linear function of the Green-Lagrange elastic strain according to S = LE (10) where L stands for the fourth order elastic sti ness tensor of the single crystal. The crystallographic nature of the plastic deformation of single crystals is described by two orthogonal unit vectors which de ne the slip system k, s and m , that stand for the slip direction and slip plane normal, respectively and that remain invariant in the intermediate con guration (Fig. 3). The available slip systems of a given crystalline material, k = 1; 2; :::; n are deter- mined by its lattice. For instance, FCC materials presented 12 slip systems characterized by the f111g planes and the <110> directions. Rice (1971) proposed the conventional ow rule of a single crystal based on this geometrical description according to X X p k k k k k L = _ s m = _ Z (11) k k k k where are internal variables of the model ( 2 q) that account for the k k accumulated plastic slip in each slip system k. The dyadic product s is the non-symmetric Schmid tensor of the system k, Z . The evolution of is normally assumed to be dependent on the internal variables chosen to describe the state q and on the stress through the resolved shear stress k k k _ = _ ( ; q): (12) The resolved shear stress is de ned in the framework of nite strains as (Cuitino ~ and Ortiz, 1992) k eT e k k = F F S : s m (13) e e which can be simpli ed in the usual case of small elastic strains (F F I) as 15 k k k k S : s m = S : Z (14) The last two ingredients of a CP model are the particular functions that dictate the shear rate _ of each slip system (eq. (12)) as well as the evolution of the internal variables q during deformation. Many di erent avours of CP have been developed in the last decades depending on these functions and the most relevant ones will be reviewed below grouped in three categories: phenomenological, physically-based and strain gradient plasticity models. 3.1. Phenomenological crystal plasticity models Crystal plasticity, in contrast with classical macroscopic plasticity, has a clear physical basis and always includes explicit microscopic information of the material, such as the geometrical de nition of the active slip systems. However, CP models can use phenomenological expressions to de ne both the slip rates and the evolution of the internal variables. This phenomeno- logical approach is based on the classical constitutive equation theory, and the internal variables that determine the state of the crystal, as well as their evolution laws, are not directly related to the microscopic physical magni- tudes or processes. The rst approaches (Rice, 1971; Hill and Rice, 1972) were rate-independent but they were relatively complex and prone to con- vergence problems because multiple combinations of shear increments can lead to the same plastic strain (Peirce et al., 1982b). Although most of the problems associated with rate-independent formulations have been overcome (Anand and Kothari, 1996), viscoplastic formulations became a very popular alternative since they were introduced by Peirce et al. (1983). The set of internal variables, q, contains the information about the accu- k k mulated plastic slip, , and the critical resolved shear stress (CRSS), g , in each slip system. For a given state, g corresponds to the minimum value of the resolved shear stress, eq. (13), to activate the plastic ow in that system. The slip rate, eq. (12), is given by a power-law function according to 1=m j j k k _ = _ sign( ) (15) where _ and m stand for the reference strain rate and the strain rate sensi- tivity parameter, respectively. The strain rate sensitivity parameter is equiv- alent to 1=n, where n stands for the strain rate sensitivity exponent. Note that only positive values of _ were considered in the original work of Peirce 16 k et al. (1983) by including positive and negative slip directions s while in eq. (15), as in most of the actual formulations, _ can take both signs and only one direction of s is considered. Regarding strain hardening, the model by Peirce et al. (1983) was de- veloped for monotonic loading and therefore only considered isotropic hard- ening. The CRSS on each system included the explicit contributions of the slip in the same system (self hardening) and of the slip on all the other slip systems (latent hardening) and the evolution of the CRSS can be written as g = h j _ j (16) kj where h are the latent hardening moduli, with h the self hardening mod- kj kk ulus. The model proposed by Peirce et al. (1983) de ned the hardening modulus, eq. (16), as h = h( )[q + (1 q) ] (17) kj kj where q is a parameter de ning the latent hardening, is the accumulated shear on all the slip systems = j _ jdt (18) and the function h( ) proposed for FCC single crystals takes the form h( ) = h sech (19) s 0 where h is the initial hardening modulus and and stand for the initial 0 0 s and the saturation values of CRSS, respectively. Alternative expressions of h( ) have been proposed to account for non- vanishing hardening rates at large plastic strains. This is the case of the Voce hardening law proposed in Tom e et al. (1984) that de nes h as h h h 0 0 0 h( ) = h 1 exp + ( + h ) exp (20) s s 0 s s 0 s 0 0 where a new parameter h is introduced to de ne the hardening slope at large plastic strains. An alternative hardening model, developed by Bassani 17 and Wu (1991), proposed a similar non-vanishing hardening rate through a parameter h using eq. (19) as starting point. The corresponding hardening moduli can be expressed as 2 h h j 0 s h = (h h )sech + h 1 + f tanh kk 0 s k s kj j6=k s 0 k (21) h = qh kj kk where h and and h stand for the self hardening and the latent hardening kk kj coecients in eq. (16). In contrast to Peirce et al. (1983), the dependence of the self hardening coecient on the slip accumulated in the di erent slip systems is not uniform in this formulation, but given by some interaction coecients f that depend on the nature of the dislocation junctions between kj the slip systems. Most of the phenomenological CP models developed since then for monotonic loading are based on either of the three hardening models summarized above. The extension of this phenomenological framework to account for cyclic deformation is done by the introduction in the constitutive model of the e ect of a backstress to formulate kinematic hardening laws at the crystal level (M eric et al., 1991; Cailletaud, 1992; Hu et al., 1992) and similar approaches have been used to model creep (Hasija et al., 2003; Venkatramani et al., 2007). The plastic shear rate depends in this case on both the CRSS and the backstress according to (Cailletaud, 1992) 1=m k k k j j g k k k k k _ sign( ) if j j g _ = (22) k k k 0 if j j < g where stands for the backstress of system k and K is a numerical param- eter. The kinematic hardening is de ned by the evolution of the backstress and Cailletaud (1992) proposed the following expression for FCC metals k k k _ = c _ dj _ j (23) where c and d are material parameters that de ne the hardening rate. Other expressions of non-linear phenomenological kinematic hardening laws at the crystal level have been proposed based in most cases in well- stablished relations at the macroscale. One of the most common expressions 18 is the Frederick-Armstrong law, adapted from macroscopic plasticity to the crystal level in M eric et al. (1991) according to k k k k _ = c _ dj _ j : (24) to account for the Bauschinger e ect in single crystal Ni-based superalloys. More recently, Cruzado et al. (2017) presented a phenomenological CP model for cyclic loading that includes the e ect of cyclic softening and an alternative evolution of the backstress to account for the mean stress relaxation. The plastic slip rate in this model is expressed as 1=m k k j j k k k _ = _ sign( ): (25) The evolution of the backstress is obtained as a simpli ed version of the Ohno-Wang macroscopic model (Ohno and Wang, 1993) limited to the rst two terms j j k k k k _ = c _ d j _ j (26) c=d where c and d are the parameters of the Frederick-Armstrong model while r is an extra parameter that controls the mean stress relaxation velocity. Finally, the cyclic softening was accounted for through a new internal variable, the cyclic plastic slip, , is given by, cyc Z Z t t X X k k = j _ jdt _ dt : (27) cyc 0 0 k k k k which was taken into account in the evolution of the CRSS, g _ = g _ ( ), cyc using a Voce type law with negative slope. 3.2. Physically-based crystal plasticity models Physically-based CP models (as opposed to phenomenological ones) con- tain a stronger physical connection with the microscopic mechanisms of plas- tic deformation. Thus, some microscopic physical quantities are included as internal variables (i.e. dislocation densities) and/or rate equations are based on the active microscopic deformation mechanisms. In addition, these mod- els usually include as input some additional microstructure information such 19 as the initial dislocation density, the volume fraction of precipitates or second phases and, very often, include the e ect of temperature. The rst ingredient of physically-based models is the relationship be- tween the plastic slip rate and the dislocation movement. This is introduced through the Orowan equation (Orowan, 1934) that connects the plastic slip k k rate on a given slip system, _ , with the mobile dislocation density, , the k k Burgers vector, b , and the average dislocation velocity, v , according to k k k k _ = b v : (28) This equation replaces the power law rate equations in the phenomeno- logical rate-dependent models, i.e. eqs. (15), (22), (25). The driving force for the dislocation movement in eq. (28) is the resolved shear stress, eq. (13), that is introduced through the average velocity. In metals with compact lat- tices (FCC and for some slip systems in HCP), lattice friction is negligible and slip occurs at very low stresses. In this case, a linear viscous relation can be postulated between the resolved stress and the velocity of a single dislocation (Hirth and Lothe, 1982) leading to k k v = (29) where B is the drag coecient, a material parameter that depends on tem- perature. If the density of the dislocations is very low, the average dislocation velocity in eq. (28) can be obtained directly from the CRSS using eq. (29), and a linear viscous relationship is established between the CRSS and the plastic slip rate. However, this relationship is not realistic in most cases be- cause dislocations have to overcome di erent barriers during slip, leading to an average dislocation velocity di erent from the one given by eq. (29). The barriers to dislocation movement can be classi ed as temperature- dependent (thermal) or independent (athermal), depending on whether ther- mal activation can help to overcome the obstacle. For instance, long-range elastic interactions among dislocations introduce an athermal threshold CRSS for the dislocation movement . Taylor (1938) determined this threshold stress for pure metals, which is given by / b (30) where is the shear modulus and 1= stands for the average distance between dislocations. 20 The thermal barriers are due to the short range interactions of dislocations with other dislocations (jogs created by the intersection of forest dislocations and their movement by vacancy generation) and point defects. The strength of the barrier at 0 K is given by and the CRSS necessary to overcome the barrier is given by + . At nite temperatures, thermal energy helps the a t dislocation to jump over the barrier, and the average dislocation velocity (eq. 28) becomes dependent on the temperature. The in uence of temperature on dislocation slip under short range inter- actions was studied in detail in Kocks (1975) with the framework of the tran- sition state theory and this work is the basis of most temperature-dependent physically-based CP models in the literature (Kothari and Anand, 1998; Ma and Roters, 2004; Cheong and Busso, 2004; Rodr guez-Gal an et al., 2015; Shahba and Ghosh, 2016). In summary, the average dislocation velocity to be inserted in eq. (28) can be written as 0 if 0 j j v = (31) k k l exp if < j j < + 0 a a t kT where l is the average distance between the obstacles in the slip system k, the attempt frequency, k the Boltzmann constant and T the absolute temperature. G( ) stands for the Gibbs free energy that has to be supplied by thermal uctuations to overcome the obstacle, which depends on the applied shear stress , and the exponential term expresses the probability of the occurrence of a jump over a short-range barrier. The evolution of the Gibbs free energy with the applied shear stress in the presence of a general array of obstacles in the slip plane can be expressed as (Kocks, 1975), j j G( ) = F 1 (32) where F is the activation free energy necessary to overcome the obstacles without the aid of an applied shear stress, and < x > stand for the Macaulay brackets, which return x if x > 0 and 0 otherwise. Finally, p and q (in the range 0 p 1 and 1 q 2) are two parameters that de ne the strength of the obstacle as a function of the distance propagated by the dislocation. A simpli cation of eq. (32) often found in CP models (Ma and Roters, 2004; Dunne et al., 2007) assumes that the obstacle strength is constant (p = q = 1) and eq. (32) can be written as 21 k G( ) = F V (33) where V is the activation volume, a material constant that determines the actual volume of the material a ected by the short range dislocation-obstacle interaction. It should be nally noted that the most important contribution to the thermal strength, , in metals with non-compact lattices (such as BCC or low density planes in HCP crystals) is the lattice friction. In this case, at 0 K is the Peierls stress, that can be obtained from the simple model by Peierls and Nabarro (Peierls, 1940; Nabarro, 1947) or from atomistic or ab initio simulations (see, for instance, Stukowski et al. (2015) for BCC W and Yasi et al. (2009) for pyramidal slip in HCP Mg). In addition to a physically-based model of _ , many physically-based CP models incorporate micromechanical internal variables, such as dislocation densities, to account directly for the strain hardening from physical consid- erations. The relationship between dislocation densities and the athermal strength is based on the Taylor model, eq. (30). This expression is usu- ally enriched by including several terms that account for the interactions between dislocations to obtain the ow stress in the di erent slip systems (here denoted as g to emphasize the relation with the CRSS introduced in the phenomenological models) . In the case of FCC crystals, this can be expressed as (Franciosi and Zaoui, 1982) k k g = b a (34) kj where is dislocation density in the slip system j and a are a set of kj non-dimensional coecients that determine the self-hardening and the la- tent hardening due to the interactions among dislocations (similar to the phenomenological expression, eq. (16)). The values a for di erent crys- kj tal lattices can be obtained from dislocation dynamics simulations (Devincre et al., 2008; Bertin et al., 2014). The dislocation densities in eq. (34), , are introduced as local internal variables and correspond to the average of dislocation length per unit volume at each point of the crystal. The evolution of these internal variables is nor- mally based in the Kocks-Mecking model (Mecking and Kocks, 1981; Estrin and Mecking, 1984). This model considers that hardening is controlled by the competition between storage and annihilation of dislocations and that 22 both processes are additive. Thus, the corresponding evolution law for one slip system k is given by _k _k = 2y j j: (35) `( ) where the rst term, 1=`( ), is athermal and controls the storage of dis- locations and `( ) is the dislocation mean free path, which corresponds to the distance travelled by a dislocation segment before it is stopped by an obstacle. In the absence of precipitates or other obstacles to the dislocation motion, `( ) = k = , where k is a material constant. The second term, 1 1 2y , is associated with the dislocation annihilation due to dynamic recov- ery. It depends on the temperature, and is characterised by y , a constant that depends on the critical annihilation distance between dislocations. It should be noted that eq. (35) can also be extended to alloys with a distribu- tion of impenetrable obstacles. In this case, the dislocation mean free path will be determined by the obstacle spacing (Estrin and Mecking, 1984). The generalization of the Kocks-Mecking model for multiple slip systems reads k k _k = 2y j _ j (36) k 1 n ` ( ; ; ) and ` q : (37) j6=k Several modi cations of the Kocks-Mecking model can be found in the lit- erature based the results obtained from simpli ed dislocation mechanics mod- els (Cheong and Busso, 2004) or dislocation dynamics simulations (Devincre et al., 2008; de Sansal et al., 2010). Moreover, the introduction of the dis- tance to grain boundary as an upper bound to the dislocation mean free path has been successfully used to simulate the e ect of grain size on the strength of Cu polycrystals (Hauoala et al., 2018). 3.3. Strain gradient crystal plasticity models The CP framework presented above is local, i.e. the material response at a given point depends only on the local values of both state and inter- nal variables at that point. Therefore, the constitutive equations are size- 23 independent because there are not length scales involved. However, the ex- perimental evidence as well as dislocation dynamics simulations show that the strength of single crystals is size-dependent when they are subjected to homogeneous (El-Awady, 2015) and inhomogeneous plastic deformation, such as nanoindentation (Stelmashenko et al., 1993; S anchez-Mart n et al., 2014a) or bending tests of single crystal cantilever beams (Motz et al., 2005; Kiener et al., 2008; Gong and Wilkinson, 2011). The development of plastic strain heterogeneities in the crystal deformation is also the reason behind the well-known dependency of the plastic strength of polycrystals on grain size (Hall, 1951; Petch, 1953). The e ect of the plastic heterogeneities in the mechanical response of a crystal was rst rationalized by Nye (1953) and Ashby (1970) as a result of the interaction between statistically stored dislocations (SSDs), which evolve from random trapping processes during plastic deformation, and geometri- cally necessary dislocations (GNDs) induced by the presence of plastic strain gradients. From the modelling viewpoint, the most common approach to ac- count for the e ect of plastic strain heterogeneities in the material response is by introducing the in uence of some plastic strain gradient measure in the constitutive equation, leading to the so-called strain gradient plasticity. This idea was introduced by Aifantis (1987) for macroscopic plasticity and further developed to account for size e ects in polycrystals (Fleck et al., 1994; Nix and Gao, 1998). This modelling framework was extended to CP by Acharya and Bassani (1995) and Shu and Fleck (1999), who in these seminal papers introduced the two di erent strategies followed in the strain gradient crystal plasticity (SGCP) models developed since then. In the rst modelling approach, lower-order SGCP, strain gradients enter only in the instantaneous hardening moduli, while the thermodynamic con- sistency is preserved (Acharya et al., 2003). This approximation allows the use of the classical mathematical framework of boundary value problems in standard plasticity (Acharya and Bassani, 1995; Dai and Parks, 1997; Busso et al., 2000; Han et al., 2005a,b; Ma et al., 2006; Cheong et al., 2005; Dunne et al., 2007). The second approach stands for the higher order SGCP models in which some internal variables are chosen as kinematic variables. This im- plies the introduction of stresses conjugated to these kinematic variables as well as the corresponding "higher order" boundary conditions (Gurtin, 2002; Gurtin and Needleman, 2005; Gurtin, 2008; Bardella, 2006; Bardella et al., 2013; Ye mov et al., 2004; Evers et al., 2004; Bayley et al., 2006; Borg et al., 2008; Niordson and Kysar, 2014). 24 The obvious bene t of lower-order formulations is their simple structure, so they can be easily implemented in existing general-purpose nite element codes. In addition, they avoid the additional higher-order boundary condi- tions, which are not always easy to interpret physically. Due to these facts, lower-order SGCP models are the basis of most physically-based CP models that account for size e ects. In addition, they have been coupled to vacancy di usion models to account for the e ect of dislocation climb, a critical pro- cess to simulate creep deformation (Geers et al., 2014). On the other hand, the main limitation of lower-order formulations is the impossibility of ac- counting for boundary layers in constrained plastic ow, as the development of these layers requires the use of higher order boundary conditions on the plastic slip elds. In this sense, higher order models are capable of account- ing for the gradient development due to the presence of a passivation layer or a grain boundary within a polycrystal. 3.3.1. Lower-order SGCP The starting point in lower-order SGCP models is the de nition of a mea- sure for the plastic strain gradients. This is usually based on the relationship between the plastic slip gradients and the density of GNDs. For single slip, the plastic slip gradient gives rise to the development of a density of GNDs, , to maintain continuity in the crystal according to (Nye, 1953) GND 1 @ = (38) GND b @x where is the plastic slip, b the Burgers vector and x stands for the coordi- nate normal to the slip direction. The generalization of the previous equation to multiple slip systems is expressed through the Nye's tensor, , which was introduced by Nye (1953) and generalized by Arsenlis and Parks (1999) as a a a = b t (39) ij GND i j a a where a stands for a straight dislocation segment of length l parallel to t a a with Burgers vector b , and stands for the length of segment a per GND unit volume. Nye's tensor accounts for the lattice curvature and, therefore, can be expressed as function of the plastic slip gradients in the case of small strains according to (Arsenlis and Parks, 1999) k k k = r s m (40) 25 where r represents the curl operator. In the case of nite strains, the Nye tensor is usually de ned as (Busso et al., 2000; Ma et al., 2006) = r F (41) although other tensorial measures of the plastic incompatibility are de ned in the literature, as the nite strain geometric dislocation tensor G proposed in Cermelli and Gurtin (2001) and given by p p G = F (r F ) (42) In lower-order gradient models, the plastic strain gradient is included in the hardening expression. This is normally done adding to the Taylor hard- ening model (eq. (34)) that relates the CRSS with the dislocation density, a term that accounts for the GND density and depends on some plastic gradient measure. Along these lines, a rst model inspired on physical considerations, was presented by Han et al. (2005a,b), who de ned a GND density for each slip system, , based in a projection of the dislocation tensor (eq. 42) in the system. Then, the e ective gradient-dependent CRSS, g , was given by eff k k g = g (g =g ) + l (43) 0 0 eff where g is the CRSS of the system k in the absence of gradients and l stands for an intrinsic length scale. It should be noted that eq. (43) resembles the phenomenological approaches of macroscopic gradient plasticity developed by Nix and Gao (1998) but the physical origin of l is clearer and can be obtained from physical considerations. If the density of GNDs in each system is explicitly resolved, the CRSS can be de ned as (Cheong et al., 2005) j j SSD GND g / b a + a (44) kj SSD kj GND j j where stands for the SSDs density and corresponds to the GNDs SSD GND SSD GND density in system j. The coecients a and a de ne the latent harden- kj kj ing interactions among the di erent slip systems. The distribution of GND densities on the di erent systems is not uniquely determined by the SSD dislocation tensor and additional constraints have to be imposed. Arsenlis and Parks (1999) proposed to minimize either the dislocation density or the 26 dislocation length to obtain the actual GND distribution. Alternatively, the GND evolution can be determined by integrating some evolution laws ob- tained by expressing the Nye tensor in terms of the spatial gradient of the slip rate (Busso et al., 2000; Ma et al., 2006; Cheong et al., 2005; Dunne et al., 2007). In order to implement these models in a nite element framework, the slip gradients or Nye tensor at the integration point level have to be determined. This task has been traditionally done using a very ecient local element approach (Busso et al., 2000; Dunne et al., 2007) in which the internal vari- ables included in the gradient term are extrapolated from the integration points to the nodes in each element. The gradients are obtained by deriv- ing the element shape functions. However, it has been recently shown by Rodr guez-Gal an et al. (2017) that this simple approach might present con- vergence problems, and alternative methods to compute the gradients based on recovery techniques are more reliable (Han et al., 2007; Rodr guez-Gal an et al., 2017). 3.3.2. Higher-order SGCP There are several alternative formulations of higher-order SGCP models. One of the most extended ones was proposed by Gurtin (2002); Gurtin and Needleman (2005); Gurtin (2008). The principal concepts and equations of this theory will be reviewed here because they are common to many other gradient models (Bardella, 2006; Bardella et al., 2013; Borg et al., 2008; Niordson and Kysar, 2014). In this theory, the plastic slip in each system, , is included in the constitutive equation as a kinematic variable (independent state variable). A central point to the model is that the work associated with each independent kinematic process should be accounted for in the energy balance. Therefore, microforces conjugated with the slip and slip gradients appear in the formulation. In the absence of body forces, the virtual power principle for a domain with boundary @ under external macroscopic surface traction t and a microscopic surface traction for each slip system reads Z Z Z Z X X k k k k k k ~ ~ ~ tv ~dA+ _ dA = : grad(v ~)dV + ( _ + grad( _ ))dV (45) where the left- and right-hand sides of this equation correspond to the ex- ternal and the internal power, respectively, and the elds v ~ and _ are the 27 k k virtual elds of the velocity and slip rate, respectively. and are the higher order stresses conjugated with the slip rate (scalar microforce) and the gradient of the slip rate (vector microforce). From the virtual power ex- pression, the resulting balance equations are the classical macroscopic force and momentum balances and the microforce balance, k k k div + = 0 for k = 1; ::; n: (46) Note that the macroscopic stress tensor enters in the microforce balance k k through the resolved shear stress , eq. (13). In this framework, can be viewed as an internal resistance force to the slip caused by the other disloca- tions and the microforce vector represents the interaction of dislocations through surfaces. The second ingredient of the theory is the constitutive equation for the microstresses. In Gurtin (2002), the elastic free energy, eq. (6), is augmented by a defect energy that depends on the dislocation tensor G eq. (42), de ned in Cermelli and Gurtin (2001). If a quadratic expression is chosen, the defect energy is given by = jGj (47) where is an scalar material parameter, with force dimensions, that can be split in a reference modulus of strength with dimensions of stress and the square of a characteristic length scale l (Bardella, 2006), = l : (48) From the expression of the defect energy, eq. (47), a linear relation is found between the microforce vector and the dislocation tensor, k 1 e k k = J F (m Gs ): (49) Finally, higher order boundary conditions should be applied in the ex- ternal boundaries. Two type of microscopic boundary conditions are usually considered in higher order SGCP, microfree and microhard boundary condi- tions. A surface S has microfree boundary conditions if n = 0 on S; k = 1; 2; :::; n (50) and this implies that slip through that surface S is not restricted. Thus, the microforce vector is directly linked to the applied macroscopic traction 28 on that boundary. The microhard boundary condition establishes that the plastic slip on a surface S is restricted, _ = 0 on S; k = 1; 2; :::; n (51) This condition emulates, e.g. a passivation layer that does not allow the dis- location ux. Other special boundary conditions can be applied at the grain boundaries based on the ow of Burgers vector at and across the boundary surface, as discussed by Gurtin and Needleman (2005). The nite element implementation of higher-order SGCP models is a complicated task. Firstly, it requires higher order continuity associated with higher order strain gradient terms, and this condition is normally addressed by means of mixed nite elements formulations (Shu and Fleck, 1999). Sec- ondly, the application of standard implicit integration algorithms to higher- order theories presents diculties in both eciency and accuracy for some relevant boundary value problems, as shown in the implementation of Gurtin macroscopic theory of 2005 proposed by Lele and Anand (2008). Di erent numerical alternatives have been used to improve the eciency of nite ele- ment implementations, using discontinuous Galerkin (Ostien and Garikipati, 2008), explicit integration (Borg et al., 2008), or the recent implicit vis- coplastic approach proposed by Panteghihi and Bardella (2016) using a spe- cial viscoplastic potential. In addition, the implementation of higher-order SGCP models implies the additional computational cost of introducing a very large number of kinematic variables. For instance, the number of kinematic variables per node will increase from 3 to 15 if a theory similar to the one presented by Gurtin (2002) is implemented in a three dimensional model of a FCC material. Thus, higher-order SGCP models have not been applied to simulate the behavior of complex three dimensional RVEs of polycrystals and most of the computational homogenization studies in the literature using SGCP are based in lower-order theories. An interesting and promising alter- native to overcome some of the limitations of the higher-order SGCP models is the use of an FFT-based framework to solve the boundary value problem of higher-order models, as shown by Lebensohn and Needleman (2016), who used the small strain version of Gurtin (2002) theory to analyze grain size e ects in polycrystals. 3.4. Deformation by twinning In addition to dislocation slip, plastic deformation in some metals with low symmetry crystal structures (such as HCP Ti, Mg and Zr) can occur 29 by twinning (Mahajan and Williams, 1973). This deformation mechanism is present when there are not enough slip systems for an arbitrary shape change of the crystal and twinning provides an additional mechanism to accommodate deformation. A mechanical twin formally corresponds to a sheared volume for which the lattice orientation is transformed into its mirror image across a so-called twin or habit plane (oblique dividing plane de ned by the twinning direction). The sheared region of the crystal undergoes an irreversible shear deformation whose value is determined by the lattice geometry and twin plane. Mechanical twinning is a process that involves two steps. The rst one is the nucleation and propagation of a thin twin band across the grain, starting normally from a grain boundary. Afterwards, the twinned region propagates in the direction perpendicular to the twin plane and eventually occupies most of the parent grain (Mahajan and Williams, 1973). The geometrical description of deformation twinning was stablished in 1965 (Bilby and Crocker, 1965). Since then, many formulations have been published, in which the geometrical description and the driving forces for twinning have been studied for di erent materials (Christian and Mahajan, 1995). The introduction of deformation by twinning in CP models was taken into account since the early developments of CP (Houtte, 1978; Tom e et al., 1991), specially due to its relevance on the deformation of HCP metallic alloys. The mechanical process of twinning is very complex and strong sim- pli cations are needed for the introduction in a CP framework. The rst models were applied in the context on mean- eld approximations (Houtte, 1978; Tom e et al., 1991) and one of the main challenges of introducing twin- ning was accounting for the large number of crystal orientations that appear due the formation of twins within each grain. Houtte (1978) tracked the evolution of the volume fraction of the twinned regions in each grain, but re- orientation was performed only on selected full grains following a statistical criterion based on the evolution of volume fraction of the twinned regions in the grain and in the entire polycrystal. Thus, the number of grains was kept constant in this approach. An improvement with respect to this approach was proposed by Tom e et al. (1991) leading to the predominant twinning reorientation (PTR) model that has been extensively used within the frame- work of the VPSC formulation to account for deformation twinning (Proust et al., 2007). In addition, an alternative model was also proposed in Tom e et al. (1991) based on a discretization of the orientation space SO3 and the representation of the orientations as volume fractions in the discrete space. 30 With this framework, the accommodation of deformation by twinning does not increase the number of orientations in the model and non-predominant orientations can also be represented. Nevertheless, the main limitaiton of these models is that they were developed within the context of mean- eld approximations for polycrystal plasticity and, therefore, cannot be easily in- cluded in single crystal plasticity formulations. The rst model that was able to describe twinning deformation of the individual single crystals was developed by Kalidindi (1998), and other sim- ilar models have been developed afterwards (Staroselsky and Anand, 2003; Kowalczyk-Gajewska, 2010; Abdolvand and Daymond, 2013a,b; Chang and Kochmann, 2015; Mareau and Daymond, 2016). Kalidindi's model has been extensively used (Zhang and Joshi, 2012; Herrera-Solaz et al., 2014b,a; Hidalgo- Manrique et al., 2015; Khan et al., 2016; J.Jung et al., 2017), and will be brie y described here. The twinning model of Kalidindi (1998) is a two-scale model. Each mate- rial point has a substructure and is divided into two phases, a parent region and a twinned region (Fig. 3), which might be potentially formed by N tw subregions. Each subregion belongs to a given twinning system and its volume fraction is f . Thus, the parent region volume fraction is given by tw 1 f . The material point is considered a composite material in which =1 the iso-strain hypothesis holds (F and F are the same in all regions). The plastic deformation is the result of three mechanisms, and the plastic ve- locity gradient in the intermediate con guration contains three terms, the standard one for dislocation slip, eq. (11), and two new contributions. The rst extra contribution, L , stands for the rate of deformation due to the tw twin transformation of a volume fraction df of the parent phase tw L = f s m (52) tw tw tw tw =1 where f = df =dt is the transformation rate in the twin system , m and tw s are unit vectors along the twin plane normal and the twining direction, tw respectively, and stands for the characteristic shear strain of the twin- tw ning mode. It is interesting to note that, contrary to dislocation slip, the accumulated plastic deformation by twinning is limited, and the maximum plastic shear that can be by accommodated is that corresponds to the tw full transformation of a material point into a given twin variant. The second contribution to plastic slip corresponds to the slip in the 31 Figure 3: Multiplicative decomposition indicating material point subdivision in parent and twin phases. Reprinted from Herrera-Solaz et al. (2014b) N slip systems of the transformed regions, L , which can be expressed sl tw sl tr as, N sl tw tw X X k k k L = f _ s m (53) sl tr sl sl =1 k =1 k k where s and m stand for the unit vectors in the slip and normal directions sl sl to the slip plane. They can be computed from the orientation of the slip planes k in the parent grain by means of the matrix, Q , Q = 2m m I: (54) tw tw that takes into account the rotation of the crystal due to twinning in the plane. An evolution equation de ning f , eq. (52), has to be speci ed to com- plete the ow rule for twinning. It is well accepted that the driving force for twin growth is the resolved shear stress on the twinning system. Accord- ingly, a viscous law depending on the resolved shear stress was proposed in Kalidindi (1998), which is equivalent to the expression used for slip in Asaro and Needleman (1985). Thus, h i if 0 _ _ f = f with h i = (55) 0 if < 0 where f is a reference transformation rate, m the strain rate sensitivity parameter and < x > stand for the Macaulay brackets that are introduced to account for the polar nature of twinning deformation. is the resolved stress in the twinning system and its value is given by parent = S : s m (56) tw tw parent where S denotes the value of the second Piola-Kircho stress in the parent region. Note that the stresses in the parent and twinned region are parent di erent due to the isostrain approach. In particular, S is obtained directly from the linear relation between the elastic strain in the intermediate con guration, E (common for both parent and twin phases) parent e S = L : E (57) where L stands for the fourth-rank elastic sti ness tensor of the crystal in its original orientation. The Piola-Kirchho stress tensor for the full integration point (containing the parent and all the twinned phases) in the intermediate con guration, S, is obtained in this case from the volume-averaged stress tensors in the di erent phases N N tw tw X X parent S = 1 f S + f S (58) =1 =1 where S is the stress in each of the twinned regions, and is obtained from an expression equivalent to eq. (57). 4. Viscoplastic self-consistent homogenization of polycrystals The computation of the e ective mechanical response and texture evolu- tion of polycrystalline materials using homogenization theory has a long tra- dition (Sachs, 1928; Taylor, 1938) and self-consistent approximations have been extensively used to deal with this problem. The 1-site VPSC the- ory of polycrystal deformation can be traced back to the work of Molinari et al. (1987a), who established a homogenization procedure based on an it- erative method involving the computation of integrals in ellipsoidal domains of the in nite-medium Green's function, customarily used in the solution of the PDEs governing the micromechanical response of heterogeneous materi- als. This formulation was implemented numerically by Lebensohn and Tom e (1993) taking into account the polycrystal anisotropy, leading to the rst 33 version of the VPSC code. Since its inception, the VPSC code has expe- rienced several improvements and extensions, e.g. recrystallization (Wenk et al., 1997); 2-site approximation for 2-phase polycrystals (Lebensohn and Canova, 1997); VPSC modelling of lamellar structures (Lebensohn et al., 1998); relative directional compliance interaction (Tom e, 1999); second-order linearization (Lebensohn et al., 2007); improved VPSC modelling of twin- ning using the PTR approach (Proust et al., 2007); dislocation density-based hardening models (Beyerlein and Tom e, 2008); climb and glide VPSC model (Lebensohn et al., 2010); dilatational VPSC for porous polycrystals (Leben- sohn et al., 2011); lattice rotation rate uctuation calculation (Lebensohn et al., 2016); improved VPSC for arbitrarily low rate sensitivities (Knezevic et al., 2016b); improved hardening laws for strain-path changes (Wen et al., 2016); VPSC prediction of intragranular misorientation evolution (Zecevic et al., 2017b), etc. The VPSC homogenization strategy is nowadays exten- sively used to simulate plastic deformation of polycrystalline aggregates and for interpretation of experimental results in metals, minerals and polymers. The self-consistent theory is one of the most commonly used homogeniza- tion methods to estimate the mechanical response behavior of polycrystals and was originally proposed by Hershley (1954) for linear elastic materi- als. For nonlinear aggregates (as those formed by grains deforming in the viscoplastic regime), several self-consistent (SC) approximations have been proposed. They di er in the procedure used to linearize the non-linear me- chanical response at the grain level, but all of them eventually make use of the original linear SC theory. They include the secant (Hill, 1965b; Hutchin- son, 1976), the tangent (Molinari et al., 1987b; Lebensohn and Tom e, 1993) and the ane (Masson et al., 2000) rst-order approximations, which are based on linearization schemes that use the information on eld averages at the grain level, and disregard higher-order statistical information inside the grains. However, the above assumption may be questionable in materials which present strong directionality and/or large variations in local properties. This is the case of low rate-sensitivity materials, aggregates made of highly anisotropic grains, voided and/or multiphase polycrystals. In all these cases, strong deformation gradients are likely to develop inside grains because of di erences in properties with neighbour crystals or phases, including voids. More accurate nonlinear homogenization methods were developed, mainly due to the work of Ponte Castaneda ~ and collaborators, to overcome the above limitations. These methods use linearization schemes at grain level that also incorporate accessible information on the second moments of the 34 stress eld distributions in the grains. These more elaborate SC formulations are based on the concept of a linear comparison material, which expresses the e ective potential of the nonlinear viscoplastic polycrystal in terms of that of a linearly viscous aggregate whose properties are determined from suitably-designed variational principles. Ponte Castaneda's ~ rst variational method was originally proposed for nonlinear composites (Ponte-Castaneda, ~ 1991) and then extended to viscoplastic polycrystals (de Botton and Ponte- Castaneda, ~ 1995). It makes use of the SC approximation for linearly vis- cous polycrystals to obtain bounds and estimates for nonlinear viscoplastic polycrystals. The more recent second-order method, proposed for nonlinear composites (Ponte-Castaneda, ~ 2002), and later extended to visoplastic poly- crystals (Liu and Ponte-Castaneda, ~ 2004), uses the SC approximation for a more general class of linearly viscous polycrystals, having a non-vanishing strain-rate at zero stress, to generate even more accurate SC estimates for viscoplastic polycrystals. The implementation of a fully anisotropic second- order approach inside the VPSC code (Lebensohn et al., 2007) has been a necessary step towards improving its predictive capability for polycrystalline materials that exhibit high contrast in local properties. Unavoidably, this improved capability came at the expense of more complex and numerically demanding algorithms. In what follows, the VPSC formulation is rst pre- sented using the ane linearization scheme (Masson et al., 2000), and the second-order linearization procedure (Liu and Ponte-Castaneda, ~ 2004) is de- scribed next. The self-consistent formulation represents the polycrystal by means of weighted, ellipsoidal, statistically-representative (SR) grains. Each of these SR grains represent the average behavior of all the grains with a particu- lar crystallographic orientation and morphology but di erent local environ- ments. These SR grains should be regarded as representing the behavior of mechanical phases, i.e. all the single crystals with a given orientation (r) belong to mechanical phase (r) and are represented by the SR grain (r). Note the di erence between mechanical phases, which di er from each other only in terms of crystallographic orientation and/or morphology, and actual phases, which di er from each other in crystallographic structure and/or composition. In what follows, SR grain (r) and mechanical phase (r) will be used interchangeably. The weights represent volume fractions and are chosen to reproduce the initial texture of the material. Each representative grain is treated as an ellipsoidal viscoplastic inclusion embedded in an e ective viscoplastic medium. Both inclusion and medium are anisotropic. Plastic 35 deformation in the inclusion is accommodated by dislocation slip activated by a resolved shear stress. As a consequence of all the above assumptions, the representation of the polycrystalline aggregate under the SC model is non-space-resolved, and corresponds to an entire class of polycrystals with microstructures consistent with a given statistical distribution. In general, homogenization models for viscoplastic deformation assumes that the elastic strains are much smaller than the plastic ones and, thus, are neglected. The plastic deformation rates are constitutively related to stress in the current con guration using small strain kinematics, resulting in relations between Cauchy stress and velocity gradient (instead of the sec- ond Piola-Kircho stress and the deformation gradient). Once the velocity gradient is obtained, the evolution of microstructure and micromechanical variables can be calculated by integrating the velocity gradient eld, or local (grain) averages of the latter, in small time increments to update the current con guration of the material. 4.1. Local constitutive behavior and homogenization Asumming small-strain kinematics for the deformation rates in the cur- rent con guration, the macroscopic velocity gradient V applied to a poly- i;j crystalline aggregate can be decomposed into an average symmetric strain- _ _ rate E = (V + V ) and an average antisymmetric rotation-rate ij i;j j;i ij (V V ). The plastic component of the deformation is assumed to be i;j j;i much larger than the elastic part and, therefore, the ow is incompressible. The viscoplastic constitutive behavior at each material point x is described by means of a non-linear rate-sensitive equation as k k (x) = Z (x) _ (x) (59) k=1 k k k k k where _ is the local strain rate and Z (x) = s (x) m (x) + m (x) s (x) is the symmetric Schmid tensor associated with slip system k and _ stands for the local shear rate on slip system k. The shear rate on each system follows a power-law relation (Peirce, 1983) j (x)j k k _ (x) = _ sign (x) = g (x) k 0 jZ (x) : (x)j s k 0 = _ sign Z (x) : (x) (60) g (x) 36 k where g is the CRSS of slip system k , _ is the reference strain rate, n the rate-sensitivity exponent and _ and stand for the strain rate and the deviatoric part of the Cauchy stress. If the shear rates are known, the lattice rotation rate, ! (or plastic spin) associated with slip activity at a single crystal material point x is given by: p k k ! _ (x) = Z (x) _ (x) (61) k k k k k where Z (x) = s (x) m (x) m (x) s (x) is the antisymmetric Schmid tensor. Note that although eqs. (59) - (61) can be used to deal with crystal deformation by slip and twinning, only slip will be considered in the examples that follow in the context of both mean- eld and full- eld approaches, to avoid the additional complication of twinning reorientation. Moreover, the constitutive behavior described by Eqs. (59) - (61) does not consider other possible crystal deformation mechanisms and microstructural evolution processes, such as climb, grain-boundary sliding or recrystalliza- tion. Assuming that the following linear relations, i.e. approximations of the actual non-linear relations, eqs. (59) - (60), hold for the SR grain (r), it can be written (r) 0 0(r) _ _ (x) = M (x) + (62) (r) 0(r) where M and are the viscous compliance tensor and back-extrapolated strain rate (strain rate under zero stress) of grain (r), respectively, which depend on the linearization assumption. For example, under the ane lin- earization, they are given by n 1 k k(r) k(r) k(r) 0(r) Z jZ : j s s s (r) M = n _ (63) k(r) k(r) g g k=1 k k(r) 0(r) jZ : j 0(r) k(r) 0(r) _ = (1 n) _ sign Z : (64) k(r) k=1 where the index (r) on the magnitudes on the right-hand side of these equa- (r) (r) 0(r) 0 0 tions indicates average over the SR (r), e.g. = h (x)i = h i . Following Hill (1965b) and Hutchinson (1976), the homogenized behavior of a linear heterogeneous medium whose local behavior is described by eq. 37 (62) can also be described by an analogous linear relation at the e ective (macroscopic) level 0 0 _ _ E = M + E (65) where E and are the overall (macroscopic) deviatoric strain rate and stress tensor, respectively, and M and E stand for the viscous compliance tensor and the back-extrapolated strain rate, respectively, of an a priori unknown Homogeneous Equivalent Medium (HEM). The response of this HEM is obtained by the linear SC method. The problem underlying the SC (r) 0(r) method is that of an inhomogeneous domain (r) of moduli M and _ , embedded in an in nite medium of moduli M and E . Invoking the concept of the equivalent inclusion (Mura, 1987), the local constitutive behavior in grain (r) can be rewritten as 0 0 _ (x) = M (x) + E + _ (x) (66) where _ (x) is an eigenstrain-rate eld, which follows from replacing the inhomogeneity by an equivalent inclusion. Rearranging and subtracting eq. (65) from Eq. (66) gives ~ (x) = L _ (x) _ (x) : (67) where the symbol "" denotes local deviations from macroscopic values of the corresponding magnitudes, and L = M . Combining eq. (67) with the equilibrium condition gives: 0 m (x) = ~ (x) = ~ (x) + ~ (x) (68) ij;j ij;j ij;j ;i where and are the Cauchy stress tensor and Cauchy the mean stress, ij respectively. From the relation between the strain rate and the velocity gradient deviations, _ (x) = (v ~ (x) + v ~ (x)), and taking into account ij i;j j;i the incompressibility condition associated with plastic deformation, L v ~ (x) + ~ (x) + ' (x) = 0 ijkl k;lj ;i ij;j (69) v ~ (x) = 0 k;k where ' (x) = L _ (x) (70) ijkl ij kl 38 is a heterogeneity or polarization eld, and its divergence f (x) = ' (x) is i ij;j an arti cial external force eld applied to the material. The set of equations (69) consists of four di erential equations with four unknowns: three are the components of the deviation from the average value of the velocity vector v ~ (x), and one is the deviation of the mean stress ~ (x). A set of N linear di erential equations with N unknown functions and a polarization term can be solved using the Green function method. G (x) km and H (x) are the Green functions associated with v ~ (x) and ~ (x), re- m i spectively, which solve the auxiliary problem of a unitary volumetric force with a single non-vanishing m-component, 0 0 0 L G (x x ) + H (x x ) + (x x ) = 0 km;lj m;i im ijkl (71) G (x x ) = 0 km;k where the Green's function G (x x ) is the velocity component in the km x -direction at point x when a unit body force in the x -direction is ap- k m 0 0 plied at point x in an in nitely extended material and H (x x ) is the corresponding response in mean stress. Once the solution of eq. (71) is obtained, the velocity eld is given by the convolution integral: 0 0 0 v ~ (x) = G (x x ) f (x ) dx (72) k ki i The set of equations (71) can be solved using Fourier transforms. Ex- pressing the Green functions in terms of their inverse Fourier transforms, the di erential equations (71) can be transformed into an algebraic set ^ ^ L k G () + ikH () = j l ijkl km i m im (73) k G () = 0 k km where k and are the modulus and the unit vector associated with a point of Fourier space = k , respectively. Calling A ( ) = L , the set j l ik ijkl (73) can be expressed as a matrix product AB = C where A, B and C are given by 39 2 3 0 0 0 A A A 11 12 13 0 0 0 6 7 A A A 21 22 23 6 7 A ( ) = 0 0 0 4 5 A A A 31 32 33 1 2 3 2 3 2 3 2 2 2 ^ ^ ^ k G k G k G 1 0 0 11 12 13 6 7 2 2 2 ^ ^ ^ 6 7 k G k G k G 0 1 0 6 7 21 22 23 6 7 B = 6 7 C = (74) 2 2 2 4 5 ^ ^ ^ 0 0 1 4 k G k G k G 5 31 32 33 ^ ^ ^ 0 0 0 ikH ikH ikH 1 2 3 Using the explicit form of matrix C, 2 3 1 1 1 A A A 11 12 13 1 1 1 6 7 A A A 21 22 23 6 7 B = A C = : (75) 1 1 1 4 5 A A A 31 32 33 1 1 1 A A A 41 42 43 gives: 2 1 k G () = A ( ) (i; j = 1; 3) (76) ij ij Comparing eqs. (74) and (75), since the components of A are real functions of , so are those of k G (). This property leads to real integrals in the ij derivation that follows. Knowing the Green's function expression in Fourier space, the solution of the eigenstrain-rate problem can be obtained using convolution integrals. Taking partial derivatives to eq. (72) leads to 0 0 0 v ~ (x) = G (x x ) f (x ) dx : (77) ki;l i k;l Replacing the expression of the arti cial volumetric force eld in eq. (77), 0 0 0 recalling that @G (x x ) =@x = @G (x x ) =@x , integrating by parts, ij ij and using the divergence theorem (Mura, 1987), 0 0 0 v ~ (x) = G (x x ) ' (x ) dx : (78) ki;jl k;l ij The integral eq. (78) provides an exact implicit solution to the problem. Furthermore, it is known from Eshelby's elastic inclusion formalism that the 40 stress and strain are constant over the domain of the inclusion (r), , if the eigenstrain is uniform over an ellipsoidal domain where the sti ness tensor is uniform. This suggests the use of an a priori unknown constant polarization within the volume of the ellipsoidal inclusion and allows to average the local eld of eq. (78) over the domain and obtain an average strain rate inside the inclusion of the form Z Z (r) 0 0 (r) v ~ = G (x x ) dx dx L _ (79) ki;jl ijmn mn k;l r r (r) (r) where v ~ and _ have to be interpreted as average quantities inside the mn k;l inclusion. Expressing the Green's function in terms of the inverse Fourier transform and taking derivatives, Z Z Z (r) 2 0 0 v ~ = k G () exp [ i (x x )] ddx dx j l ki k;l r r (r) (r) L _ = T L _ : (80) ijmn ijmn mn klij mn Writing d in spherical coordinates (d = k sin dkdd') and using the relation (76), the Green's interaction tensor T can be expressed as klij Z Z T = A ( ) ( ) sin d d' (81) klij j l ki 0 0 where Z Z Z 0 0 2 ( ) = exp [ i (x x )] dx dx k dk (82) r r Integrating eq. (82) inside an ellipsoidal grain of radii (a; b; c) (Berveiller et al 1987) and replacing in eq. (81) leads to Z Z 2 1 abc A ( ) j l ki T = sin d d' (83) klij [ ( )] 0 0 1=2 2 2 2 where ( ) = [(a ) + (b ) + (c ) ] . The symmetric and antisym- 1 2 3 metric Eshelby tensors (functions of and the shape of the ellipsoidal inclusion, representing the morphology of the SR grains) are de ned as S = T + T + T + T L (84) ijkl ijmn jimn ijnm jinm mnkl P = T T + T T L (85) ijkl ijmn jimn ijnm jinm mnkl (r) (r) ~ ~ and the average strain rate and rotation rate deviations, _ and ! _ , in the ellipsoidal domain are obtained by taking symmetric and antisymmetric components to eq. (80) and using eqs. (84) - (85) (r) (r) _ = S _ (86) (r) (r) 1 (r) ~ ~ ! _ = P _ = PS _ (87) (r) (r) (r) (r) ~ _ ~ _ where _ = E _ and ! _ = ! _ are deviations of the average strain rate and rotation rate inside the inclusion, respectively, with respect to the (r) corresponding overall magnitudes, and _ is the average eigenstrain-rate in the inclusion. Therefore, the lattice rotation rate eld is given by (r) p ! _ (x) = + ! _ ! _ (x) (88) 4.2. Interaction and localization equations Taking volume averages over the domain of the inclusion on both sides of eq. (67) gives 0(r) (r) (r) ~ = L _ _ (89) and replacing the eigenstrain-rate given by eq. (86) into eq. (89), leads to the interaction equation : 0(r) (r) ~ ~ ~ _ = M (90) where the interaction tensor is given by: M = (I S) SM: (91) Replacing the constitutive relations of the inclusion and the e ective medium in the interaction equation and, after some manipulation, leads to the localization equation : 0(r) (r) 0 (r) = B + b (92) 42 where the localization tensors are de ned as: (r) (r) ~ ~ B = M + M M + M (93) (r) (r) 0 0(r) ~ _ b = M + M E _ (94) 4.3. Self-consistent equations The derivation presented in the previous sections solves the problem of an equivalent inclusion embedded in an e ective medium. This result is used in this section to construct a polycrystal model, in which each SR grain (r) stands for an inclusion embedded in an e ective medium that represents the polycrystal. The properties of such medium are not known a priori and have to be found through an iterative procedure. If the stress localization equation, eq. (92) is introduced in the local constitutive equation, eq. (62), averaged over the SR grain (r), (r) (r) (r) (r) (r) 0(r) = M B + M b + : (95) Taking volumetric average to eq. (95) and enforcing the condition that the average of the strain-rates over the aggregate has to coincide with the macroscopic quantities, leads to, (r) E = _ (96) where the brackets denote average over the SR grains weighted by the asso- ciated volume fraction. From the macroscopic constitutive relation, eq. (65), the following self-consistent equations are obtained for the HEM's compliance and back-extrapolated term: (r) (r) M = M B (97) 0 (r) (r) 0(r) E = M b + _ (98) 4.4. Linearization assumptions As indicated above, di erent choices are possible for the linearized behav- ior at the grain level, and the results of the homogenization scheme depend on this choice. Several rst-order linearization schemes are presented below, de ned in terms of the stress rst-order moment (average) inside SR grain (r). 43 The secant approximation (Hill, 1965b; Hutchinson, 1976) assumes the following linearized secant moduli: n 1 k(r) k(r) k(r) 0(r) Z jZ : j s s s (r) M = _ (99) sec k(r) k(r) g g 0(r) _ = 0 (100) sec k(r) k(r) where the index (r) in Z and g indicates uniform (average) values of these magnitudes, corresponding to a given orientation and hardening state associated with the SR grain (r). Under the ane approximation (Masson et al., 2000), the ane moduli are given by n 1 k(r) k(r) k(r) 0(r) Z jZ : j (r) s s s M = n _ (101) aff k(r) k(r) g g k(r) 0(r) jZ : j 0(r) k(r) 0(r) _ = (1 n) _ sign Z : (102) aff s k(r) In the case of the tangent approximation (Molinari et al., 1987b; Leben- sohn and Tom e, 1993), the tangent moduli are formally the same as in the (r) (r) 0(r) 0(r) ane case: M = M and _ = _ . However, Molinari et al. (1987b) tg tg aff aff used the secant compliance, eq. (99), instead of the ane one in order to avoid the iterative adjustment of the macroscopic back-extrapolated term to adjust M (to be denoted M ), in combination with the tangent-secant rela- sec tion M = nM (Hutchinson, 1976). Then, the expression of the interaction tg sec tensor is given by 1 1 M = (I S) SM = n (I S) SM (103) tg sec Qualitatively, the interaction eq. (90), indicates that the larger the in- teraction tensor, the smaller the deviation of grain stresses with respect to the average stress. As a consequence, the tangent approximation tends to a uniform stress state or lower-bound approximation for n ! 1 (Sachs, 1928). This rate-insensitive limit of the tangent formulation is an artefact created by the use of the above tangent-secant relation of the non-linear polycrystal in the self-consistent solution of the linear comparison polycrystal. On the 44 contrary, the secant interaction has been proven to tend to a uniform strain- rate state or upper-bound approximation (Taylor, 1938) in the rate insensitive limit. 4.5. Second-order formulation The more accurate second-order approximation to linearize the behav- ior of the mechanical phase provides improved micromechanical predictions, which are obtained from the calculation of average uctuations of the stress distribution inside the linearized SR grains. The methodology to obtain these uctuations was derived by Bobeth and Diener (1987),Kreher (1990) and Parton and Buryachenko (1990) for composites and extended by Leben- sohn et al. (2007) for polycrystals. The e ective stress potential U of a linearly viscous polycrystal described by eq. (68) may be written in the form (Willis, 1981; Laws, 1973), 1 1 0 0 0 0 U = M : ( ) + E : + G (104) 2 2 where G is the power under zero applied stress. The self-consistent expression for M and E , eqs. (97) and (98), can be re-written as (r) (r) (r) (r) (r) M = M B = c M B (105) X X 0 (r) (r) 0(r) (r) (r) (r) 0(r) (r) 0(r) (r) E = M b + _ = c M b + _ = c _ B (106) r r (r) where c is the volume fraction associated with SR grain (r). The corre- sponding expression for G is (r) 0(r) (r) G = c _ : b : (107) The average second-order moment of the stress eld over a SR grain (r) of the polycrystal is a fourth-rank tensor given by: 2 @U (r) T 0 0 i = (108) (r) (r) c @M Replacing eqs. (105) - (107) in (108), leads to 45 0 1 @M 1 @E 1 @G (r) 0 0 0 0 0 i = : ( ) + : + (109) (r) (r) (r) (r) (r) (r) c @M c @M c @M The algorithmic expressions to calculate the partial derivatives on the right-hand side can be found in Lebensohn et al. (2007). Once the average second-order moments of the stress eld over each SR grain (r) are obtained, the implementation of the second-order procedure follows the work of Liu and Ponte-Castaneda ~ (2004). The covariance tensor of stress uctuations in the SR grains of the linear comparison polycrystal is given by: (r) (r) 0 0 0(r) 0(r) C = h (110) and the average and the average uctuation of resolved shear stress on slip system k of SR grain (r) is given by k(r) k(r) 0(r) = Z : (111) 1=2 (r) k(r) k(r) k(r) k(r) ^ = Z : C : Z (112) s s k(r) where the positive (negative) branch should be selected if is positive (negative). The slip potential associated with slip system k of the nonlinear polycrystal is de ned as: n+1 j j k 0 ( ) = : (113) n + 1 g The linearized local behavior associated with SR grain (r) is then given by: (r) 0(r) (r) 0(r) _ = M + _ (114) SO SO with: (r) k(r) k(r) k(r) M = Z Z (115) s s SO 0 (r) k(r) k(r) = e Z (116) SO s 46 0k k where ( ) = d =d stand for two scalar magnitudes associated with each slip system k of each SR grain (r) are de ned as 0k(r) k(r) 0k(r) k(r) k(r) = (117) k(r) k(r) k(r) 0k(r) k(r) k(r) k(r) e = (118) Once the linear comparison polycrystal is de ned by eqs. (115)-(116), di erent second-order estimates of the e ective behavior of the nonlinear aggregate can be obtained. Approximating the potential of the nonlinear polycrystal in terms of the potential of the linear comparison polycrystal and a suitable measure of the error, Liu and Ponte-Castaneda ~ (2004) generated the following expression (corresponding to the so-called energy version of the second-order theory) for the e ective potential of the nonlinear polycrystal X X 0 (r) k(r) k(r) 0k(r) k(r) k(r) k(r) U ( ) = c ^ + ^ (119) from where the e ective response of the homogenized polycrystal can be 0 0 obtained as E = @U ( ) @ . The alternative constitutive equation version of the second-order theory simply consists in making use of the e ective stress-strain-rate relations for the linear comparison polycrystal. The e ective strain is obtained as X X (r) k(r) 0k(r) k(r) E = c Z (120) r k Both versions of the second order theory give slightly di erent results, depending on non-linearity and local anisotropic contrast. Such gap is rela- tively small compared with the larger di erences obtained with the di erent SC approaches. The constitutive equation version is { in principle { less rigor- ous since it does not derive from a potential function, but has the advantage that can be obtained by simply following the ane algorithm described in the previous sections, using the linearized moduli de ned by eqs. (115)-(116). 4.6. Numerical implementation To illustrate the practical use of the VPSC formulation, the steps required to predict the local and overall viscoplastic response of a polycrystal are 47 detailed below. Starting for convenience with an initial Taylor guess, i.e. (r) 0(r) _ = E for all grains, the following non-linear equation is solved to get , k(r) 0(r) Z : k(r) k(r) 0(r) E = _ Z sign Z : (121) s s k(r) and an appropriate rst-order linearization scheme is used to obtain initial 0(r) (r) values of M and _ , for each SR grain (r). Next, initial guesses for the macroscopic moduli M and E are obtained (usually as simple averages of the local moduli). The initial guess for the macroscopic stress can be obtained from them and the applied strain rate, eq. (65), while the Eshelby tensors S and P can be calculated using the macroscopic moduli and the ellipsoidal shape of the SR grains by means of the procedure described in section 4.1. Subsequently, the interaction tensor M, eq. (91), and the (r) (r) localization tensors B and b , eqs. (94) and (95), can be calculated as well. With these tensors, new estimates of M and E are obtained by solving iteratively the self-consistent eqs. (97) and (98). After achieving convergence on the macroscopic moduli (and, consequently, also on the macroscopic stress and the interaction and localization tensors), a new estimation of the average grain stresses can be obtained using the localization relation, eq. (92). If the recalculated average grain stresses are di erent from the input values within certain tolerance, a new iteration should be started, until convergence is reached. If the chosen linearization scheme is the second-order formulation, an additional loop on the linearized moduli is needed, using the improved estimates of the second-order moments of the stress in the grains, obtained by means of the methodology described in section 4.4. This additional loop roughly increases the calculation time by one order of magnitude with respect to rst-order linearizations. When the iterative procedure is completed, the average shear rates on the slip system (k) in each grain (r) are calculated as: k(r) 0(r) Z : k(r) k(r) 0(r) _ = _ sign Z : (122) k(r) These average shear rates are in turn used to calculate the lattice rotation rates associated with each SR grain: (r) (r) p(r) _ ~ ! _ = + ! _ ! _ (123) where, following eq. (61), 48 X p(r) k(r) k(r) ! _ = Z _ (124) k(r) where Z is the (uniform) antisymmetric Schmid tensor of slip system k in SR grain (r). The above numerical scheme can be used to predict the texture devel- opment by applying the viscoplastic deformation to the polycrystal in incre- mental steps. The latter is done by assuming constant rates during a time interval t (such that Et corresponds to a macroscopic strain increment of the order of a few percents) and using: a) the strain rates and rotation rates (times t) to update the shape and orientation of the SR grains, and b) the shear rates (times t) to update the critical stress of the deformation systems due to strain hardening using one the models described in section 3, after each deformation increment. Note that the above explicit update schemes relies on the fact that the orientation and hardening variables evolve slowly within the adopted time interval. Otherwise, t should be chosen smaller. 5. Full- eld homogenization of polycrystals Full- eld homogenization aims to predict the macroscopic response and microscopic eld distribution of heterogenous materials based on the nu- merical simulation of the mechanical response of an RVE of the material microstructure (B ohm, 2004). The method is computationally expensive be- cause it involves the solution of a boundary value problem that might contain a large number of degrees of freedom. Nevertheless, it can provide more ac- curate predictions of the macroscopic response than mean- eld homogeniza- tion models (such as VPSC) because it includes more accurate information of the details of the microstructure. In addition, full- eld homogenization of polycrystals can account for the e ect of microstructural details that can- not be easily included in mean- eld models such as relative grain sizes and grain shape distribution, as well as grain boundary misorientation distribu- tions. Finally, it provides very accurate information for the local values of the stress and strain elds as well as of the state variables throughout the mi- crostructure. This information is critical for predicting damage localization and failure of heterogeneous materials. Several methods have been used to simulate the mechanical response of an RVE in within the framework of full- eld homogenization. Among them, the two most common approaches are the FFT-based method (Moulinec and 49 Suquet, 1998), extended to viscoplastic polycrystals in Lebensohn (2001), and the FEM. 5.1. Homogenization using the nite element method The rst attempts to use the FEM to obtain the mechanical response of a polycrystal were carried out in the 90s by Kalidindi et al. (1992), Bronkhorst et al. (1992) and Miehe et al. (1999). These rst works were limited to 2-D and the RVEs used for the simulations were just a regular arrangement of - nite elements in which each element represented a whole grain. These studies were able to capture the macroscopic response and texture evolution under large strains but thedistribution of local elds was missing due to the rudi- mentary representation of the microstructure. Indeed, these simulations were closer to mean- eld polycrystalline models (such as Taylor or VPSC) than to full- eld simulations because the elds in each grain were constant (or just vary linearly) and did not take into account the strong gradients found in actual polycrystals. The use of complex RVEs with more realistic representa- tions of the polycrystalline microstructure started at the early 2000s (many years after the development of CP theory and of the rst numerical imple- mentations) due to the high computational cost of these type of models. The rst 3-D simulations used regular arrangements of grains (for example, rhom- bic dodecahedra in Mika and Dawson (1999)) but each grain was represented using several nite elements and, thus, the strong deformation gradients due to the accommodation of the plastic strain incompatibility between adjacent grains were accounted for. Finally, the use of several elements per grain was combined with more complex representation of the microstructure in 3- D, leading to the rst realistic full- eld simulations (Barbe et al., 2001a,b). Since then, full- eld simulation using the FEM has been a very active area of research. The nite element implementation of CP models has been reviewed by Roters et al. (2010b,a) and therefore will not presented here. Only some notes on the boundary conditions used in the boundary value problem will be reviewed due to their implication in obtaining the macroscopic elds as result of the full eld simulations. The boundary conditions used in the simulation of the RVE deforma- tion are a key issue in full- eld homogenization. These conditions are en- ergetically consistent if the stress and strain microscopic elds and the cor- responding macroscopic elds ful ll the so-called Hill-Mandel principle of macro-homogeneity (Hill, 1965b). This principle states that the macroscopic stress power should be equal to the volume average of the microscopic stress 50 power. Four di erent types of boundary conditions can be imposed to the RVE resulting in micro elds compatible with this condition (B ohm, 2004; Peric et al., 2011; Ostoja-Starzewski, 2006): (1) statically uniform boundary conditions where homogeneous surface tractions are applied to the RVE faces; (2) kinematically uniform boundary conditions where constant displacements are applied to the RVE boundary; (3) mixed conditions combining uniform tractions and constant displacements on the di erent RVE boundaries and (4) periodic boundary conditions. For any given RVE, the results obtained under periodic boundary conditions (even in the case of non-periodic microstruc- tures) lie between the corresponding lower and upper estimates obtained with homogeneous boundary conditions (1) and (2) and show a faster convergence towards the actual e ective response with increasing the RVE size (Hazanov and Huet, 1994). Moreover, periodic boundary conditions reproduce exactly the deformation of an in nite media constructed by a periodic arrangement of the RVE. For these reasons, periodic boundary conditions are the common choice for full- eld homogenization and will be described here. Periodic boundary conditions assume that the whole RVE deforms as a jigsaw puzzle. If a RVE is selected such that the periodicity of the microstruc- ture is cubic, then the periodic translation of the cell in the three directions will ll the whole space (Figure 4(a)). Let l = l e be the three orthogonal i i i vectors de ning the cubic periodicity and let e be the corresponding unit vectors de ning the basis. If x ; x ; x are the coordinates of a point in the 1 2 3 RVE in the system de ned by e , the periodic boundary conditions link the local displacement vector u of the nodes on opposite faces of the RVE with the far- eld macroscopic deformation gradient F according to, u(x ; x ; 0) u(x ; x ; L ) = (F 1)l 1 2 1 2 3 3 u(x ; 0; x ) u(x ; L ; x ) = (F 1)l (125) 1 3 1 2 3 2 u(0; x ; x ) u(L ; x ; x ) = (F 1)l : 2 3 1 2 3 1 The resulting deformed cell preserves the cubic periodicity, but the new vectors de ning the periodicity are given by Fl , as it can be observed in Figure 4(b). These boundary conditions are implemented in the FEM by coupling the displacement degrees of freedom of the nodes in the RVE surface. The most direct implementation of periodic boundary conditions requires a periodic mesh, i.e. the 2D discretization of two opposite faces of the RVE boundary 51 (b) (a) (x1,x2+l2) Fl (x1+l1,x2) ((x x1 1,x ,x2 2)) Fl e e (x1,x2) e e Figure 4: Deformation of a periodic RVE of a polycrystal under periodic boundary condi- tions in 2-D. (a) Undeformed con guration. (b) Deformed shape under biaxial deforma- tion. should be identical, although some interpolation methods were proposed in the literature to overcome this requirement (Nguyen et al., 2012). The ap- plication of the eq. (125) to the nodes of a periodic mesh can be done by constraint elimination (Peric et al., 2011; Michel et al., 1999) or by means of the Lagrange multiplier method (B ohm and Han, 2001; Segurado and LLorca, 2002). Standard nite element codes do not allow to apply direct constraint elimination and the common approach is then to use Lagrange multipliers by imposing a relation in the displacement of some nodes in the model (Segu- rado and LLorca, 2002, 2013). In particular, a master node M (i = 1; 2; 3) can be de ned for each pair of opposite cube faces and the value of the far- eld macroscopic deformation gradient is imposed to the RVE through the displacement of those master nodes according to u(M ) = (F 1)l : (126) i i If some components of the far- eld deformation gradient are not known a priori (e.g. under uniaxial tension), the corresponding e ective stresses are set instead. This task is carried out by applying nodal forces P to the 52 corresponding master node i and degree of freedom j according to P (M ) = (e ) A (127) j i i j i where A is the area of the cell perpendicular to direction e . Under small i i strains, the area A corresponds to the undeformed geometry and eq. (127) is used to obtain the value of the forces to be applied as function of the target stress. If nite deformations are applied, A corresponds to the actual area and the relation between the force and the corresponding Cauchy stress component is not known a priori because the actual value of the area changes with the deformation. In this case, to impose a given history of any stress component, an iterative procedure should be applied to correct the applied forces with the cell deformation. 5.2. Homogenization using the fast Fourier transform The large number of degrees of freedom required by full- eld homogeniza- tion limits the size of the microstructures that can be simulated using the FEM. Conceived as a very ecient alternative to the FEM, a numerical for- mulation based on the use of the FFT algorithm was originally proposed by Moulinec and Suquet (1994) to compute convolution integrals over the entire periodic RVE involving the periodic Green's function of a reference medium. In the original method, known as the basic scheme, the solution is obtained using a xed-point iterative scheme. The contrast in phase properties has a strong impact in the convergence rate (the higher the contrast, the more numerically demanding the convergence). Moreover, the properties of the reference medium, although do not in uence the nal results, also strongly a ect the convergence rate. Several modi cations to the Moulinec-Suquet scheme have been proposed afterwards (Eyer and Milton, 1999; Michel et al., 2000; Zeman et al., 2010; Monchiet and Bonnet, 2013; Kabel et al., 2014) based on modifying the iterative algorithm in order to improve the conver- gence rate. All these approaches rely on the use of a reference medium and, therefore, their performance is linked to an appropriate election of its prop- erties. An alternative approach to FFT-based homogenization in which the prob- lem is derived using a weak formulation, analogous to nite elements, was recently proposed (Vondrejc et al., 2014; de Geus et al., 2017; Zeman et al., 2017). This approach does not require a reference medium and the Green's function of the classical schemes are replaced here by projection operators 53 that ensure the compatibility of the strain elds. The similarities of this approach with Galerkin nite elements make very easy the introduction of generic non-linear constitutive equations, including their corresponding tan- gents to accelerate the convergence. FFT-based methods were originally developed Moulinec and Suquet (1994) to carry out computational homogenization in periodic composites and have been profusely used for this class of materials (Moulinec and Suquet, 1998; Eyer and Milton, 1999; Michel et al., 1999; Idiart et al., 2006; Brisard and Dormieux, 2010; Moulinec and Silva, 2014; Willot, 2015). In this type of problems, the source of heterogeneity is related to the spatial distribution of phases with di erent mechanical properties and the phase behavior is given in general by simple constitutive equations (in many cases, by isotropic linear elastic relations). The extension of FFT-based homogenization to polycrys- tals was proposed rst in Lebensohn (2001) to exploit the superior numerical performance of CP-FFT with respect to CP-FE, allowing a large reduction of the computational time and/or the use of very complex and detailed poly- crystalline RVEs. Since then, di erent numerical implementations of the FFT-based method for polycrystals have been developed, depending on the constitutive equation adopted to represent the behavior of each single crystal material point. Classic schemes have been used to solve the linearized prob- lem In almost all the cases, requiring the de nition of a reference medium. An exception is the recent work by Lucarini and Segurado (2018) who adapted the Galerkin FFT-based approach to study the cyclic behavior and fatigue performance of polycrystals. A summary of the type of problems and constitutive regimes in which FFT-based homogenization has been applied to polycrystals includes thermo- elasticity (Brenner et al., 2009), rigid-viscoplasticity (Lebensohn, 2001; Leben- sohn et al., 2005, 2008, 2009), elasto-viscoplasticity for small strains (Leben- sohn et al., 2012; Grennerat et al., 2012) and large strains (Eisenlohr et al., 2013; Lucarini and Segurado, 2018). More complex behaviors include dilata- tional plasticity (Lebensohn et al., 2013), strain-gradient plasticity (Leben- sohn and Needleman, 2016), curvature-driven plasticity (Upadhyay et al., 2016), and transformation plasticity (Richards et al., 2013; Otsuka et al., 2018). Recent implementations of eld dislocation mechanics theories (Bren- ner et al., 2014; Berbenni et al., 2014) and discrete dislocation mechanics for- mulations (Graham et al., 2016; Bertin et al., 2015) are also examples of the great potential of FFT-based formulations to provide the numerical eciency needed for implementation of these very powerful, but also very complex and 54 numerically-demanding, emerging micromechanical formulations. In what follows, the rigid-viscoplastic FFT-based formulation (VPFFT) for fully-dense, incompressible polycrystals, and the dilatational extension of the latter (DVPFFT) to predict microstructural e ects on porosity evolution in polycrystals with intregranular voids are presented in order to illustrate the methodology and provide the theory and algorithms behind one of the examples of section 7. For details of other numerical implementations of the FFT-based formulation, the reader is referred to the speci c papers listed above. Brie y, the viscoplastic FFT-based formulation consists in iteratively ad- justing a compatible strain rate eld, related with an equilibrated stress eld through a constitutive potential, such that the average of local work rates is minimized. The method is based on the fact that the local mechanical response of a heterogeneous medium can be calculated as a convolution in- tegral between Green functions associated with appropriate elds of a linear reference homogeneous medium and the actual heterogeneity eld. For a pe- riodic medium, use can be made of the discrete Fourier transform to reduce convolution integrals in real space to simple products in Fourier space. Thus, the FFT algorithm can be utilized to transform the heterogeneity eld into Fourier space and, in turn, to get the mechanical elds by transforming that product back to real space. However, the actual heterogeneity eld depends precisely on the a priori unknown mechanical elds. Therefore, an iterative scheme has to be implemented to obtain, upon convergence, a compatible strain rate eld and a stress eld in equilibrium. A periodic RVE of the polycrystal is discretized into N N N Fourier 1 2 3 points. This discretization determines a regular grid in the Cartesian space x and a corresponding grid in Fourier space . Velocities and tractions along the boundary of the unit cell are left undetermined. A velocity-gradient V (which can be decomposed into a symmetric strain rate and a antisym- i;j _ _ metric rotation rate: V = E + ) is imposed to the unit cell. The local i;j ij ij strain rate eld is a function of the local velocity eld, i.e. _ (v (x)), and can ij k _ ~ be split into its average and a uctuation term: _ (v (x)) = E + _ (v ~ (x)), ij k ij ij k where v (x) = E x +v ~ (x). By imposing periodic boundary conditions, the i ij j i velocity uctuation eld v ~ (x) is assumed to be periodic across the boundary of the unit cell, while the traction eld is antiperiodic, to meet equilibrium on the boundary between contiguous unit cells. The local constitutive relation between the strain-rate _ (x) and the devi- atoric stress (x) is given by the same rate-sensitivity relation used within 55 0 the VPSC framework, eq. (59). If a fourth rank tensor L is chosen as the sti ness of a linear reference medium (the choice of L can be quite arbitrary, but the speed of convergence of the method will depend on this choice). The polarization eld ' (x), eq. (70), is de ned as ij 0 0 ' (x) = ~ (x) L _ (x) (128) ij kl ij ijkl Then, the Cauchy stress deviation can be written as 0 m ~ (x) = L _ (x) + ' (x) + ~ (x) (129) ij kl ij ij ijkl and combining Eq. (129) with the equilibrium ( (x) = ~ (x) = 0), the ij;j ij;j incompressibility condition, and the relation _ (x) = (v ~ (x) + v ~ (x)) ij i;j j;i leads to, 0 m L v ~ (x) + ~ (x) + ' (x) = 0 k;lj ij;j ijkl ;i (130) v ~ (x) = 0 k;k This system of di erential equations is formally equivalent to system (69). However, both systems actually di er in that: a) the HEM's sti ness mod- ulus L of eq. (69) is replaced in eq. (130) by the sti ness of a linear refer- ence medium L , and b) the polarization eld in Eq. (130) has in general non-vanishing values throughout the unit cell and is periodic (owing to the RVE periodicity), while the polarization eld in Eq. (69) vanishes outside the domain of the inclusion. The auxiliary system involving periodic Green functions is, then, given by eq. (73), 0 0 0 0 L G (x x ) + H (x x ) + (x x ) = 0 km;lj m;i im ijkl (131) G (x x ) = 0 km;k After some manipulation, the velocity and the velocity gradient deviation elds are given by the convolution integrals 0 0 0 v ~ (x) = G (x x ) ' (x ) dx (132) k ki;j ij 0 0 0 v ~ (x) = G (x x ) ' (x ) dx : (133) i;j ik;jl kl Convolution integrals in direct space are simply products in Fourier space 56 ^ ^ v ~ () = ( i ) G () ' ^ () (134) k j ki ij ^ ^ v ~ () = () ' ^ () (135) i;j kl ijkl where the symbol "^" indicated a Fourier transform. The Green operator in eq. (135) is de ned as = G . ijkl ik;jl ^ ^ The tensors G () and () can be calculated by taking Fourier trans- ij ijkl form to the set of eqs. (131) ^ ^ L G () + i H () = l j km i m im ijkl (136) G () = 0 k km 0 0 De ning the 3x3 matrix A () = jL and embedding it in the 4x4 ik ijkl matrix A () as 0 0 0 A A A 11 12 13 0 0 0 A A A 21 22 23 A () = (137) 0 0 0 A A A 31 32 33 1 2 3 it leads to (see eqs. (74)-(76)) G () = A (i; j = 1; 3) (138) ij ij ^ ^ () = G () (139) j l ik ijkl Numerical implementation Assigning an initial guess values to the strain-rate eld in the regular grid d 0 d 0 d ~ _ x (e.g. _ x = 0 ) _ x = E ), and computing the correspond- ij ij ij 00 d ing stress eld (x ) from the local constitutive relation, eq. (59) (which ij k d requires to know the initial values of the critical stresses g x , and the Schmid tensors Z x ), an initial guess for the polarization eld in direct 0 d space ' x , eq. (59), can be obtained which in turn can be Fourier- ij 0 0 d transformed to obtain ' ^ . Furthermore, assuming that x = ij ij 0 d x is the initial guess for a eld of Lagrange multipliers associated with ij the compatibility constraints, the iterative procedure based on Augmented Lagrangians proposed by Michel et al. (2000) can be used. If the polarization 57 eld after iteration n is known, the n +1-th iteration starts by computing the new guess for the kinematically-admissible strain-rate deviation eld ^ d sym d d d ^ n+1 n n+1 ~ ~ d = ' ^ ; 8 6= 0; and d (0) = 0 (140) ij ijkl kl ij sym where is the Green operator, appropriately symmetrized. The corre- ijkl sponding eld in real space is, thus, obtained by application of the inverse FFT, i.e. n o 1 ^ d n+1 d n+1 ~ ~ d x = t d (141) ij ij and the new guess for the deviatoric stress eld is calculated from (omitting subindices): k d 0n+1 d Z x : x ( ) ( ) 0n+1 d 0 k d k d 0n+1 d x + L : _ Z x sign Z x : x = s s k k d g x ( ) n d 0 n+1 d _ ~ = x + L : E + d x (142) The iteration is completed with the calculation of the new guess for the Lagrange multiplier eld according to n+1 d n d 0 n+1 d n+1 d ~ ~ x = x + L : _ x d x (143) Equations (142)-(143) guarantee the convergence of _ x (i.e. the strain- rate eld related with the stress through the constitutive equation) towards d x (i.e. the kinematically-admissible strain-rate eld) to ful l compat- ibility, and of the Lagrange multiplier eld x towards the stress eld 0 d x to ful l equilibrium. Time integration Upon convergence, the microstructure can be updated using an explicit scheme. The resulting shear strain-rate eld k d 0 d Z x : x k d k d 0 d _ x = _ sign Z x : x (144) k d g (x ) 58 can be assumed to be constant during a time interval [t; t + t]. The macro- scopic and local strain increments are then calculated as E = E t and ij ij d d x = _ x t and the local crystallographic orientations are up- ij ij dated according to the following local lattice rotation d d d d ! x = x + ! _ x ! _ x t (145) ij ij ij ij d d where ! _ x can be obtained from eqs. (61) and (176), and ! _ x is ij obtained back-transforming the converged antisymmetric eld ^ antisym ^ d d d d ~ ^ ~ ! _ = ' ^ ; 8 6= 0; and ! _ (0) = 0: (146) ij ijkl kl ij The CRSSs of the deformation systems associated with each material point should also be updated after each deformation increment due to strain hardening using one the models described in section 3. After each time increment, the new position of the Fourier points can be determined calculating the velocity uctuation term v ~ x back-transforming eq. (132), and d d d d X x = x + E x + v ~ x t: (147) i ij i i j Evidently, the set of advected Fourier points no longer forms a regular grid, after the very rst deformation increment due to the heterogeneity of the medium. A rigorous way of dealing with this situation was proposed by Lahellec et al. (2001) based on the particle-in-cell method (Sulsky et al., 1995). Extension to dilatational viscoplasticity The previous FFT-based algorithm can be adapted to the case of dilata- tional viscoplasticity of voided materials with polycrystalline or homogeneous isotropic matrix. The constitutive behavior of Fourier points belonging to a void is simply = 0 . In the cases of porous materials with crystalline or isotropic matrix, the constitutive equation for Fourier points with material properties is given, respectively, by eq. (59), or by n 1 3 _ eq d d _ x = x (148) 0 0 59 where is the ow stress. Given the e ective compressibility of the material due to the presence of voids, the di erential equation whose solution is the Green's function G x associated with the velocity eld is given by eq. (131) km o 0 0 L G (x x ) + (x x ) = 0 (149) km;lj im ijkl hence (see eqs. (138)-(139)) 1 o G () = A () and A () = L (150) ij ik l j ij ijkl and ^ ^ () = G () : (151) j l ik ijkl The algorithm for a full stress tensor imposed to a unit cell representing a porous material requires an initial guess for the average strain rate 0 E 0 kk _ _0 E = E + (152) ij ij ij _ _ where E and E are the deviatoric and hydrostatic parts, which will be kk ij adjusted iteratively. Initial guess values also need to be assigned to the local strain rate eld. The deviatoric part is taken to be uniform 0 0 0 0 0 0 _ (x) = 0 ) _ (x) = E : (153) ij ij ij An initial guess of uniform dilatation is assumed for the Fourier points representing voids, o o _ (x) = 1 E (154) kk kk where f is the current porosity in the RVE. For these initial values, the corresponding stress eld in the material points (x) is obtained by in- verting the local constitutive relation given by eq. (59) (for polycrystals) or eq. (148) (for isotropic matrix), whereas for points belonging to voids, the stress simply vanishes. The initial speci cation of these elds and reference sti ness allows us to calculate the initial guess for the polarization eld in direct space, eq. (128), which in turn is Fourier-transformed, used in combi- nation with eq. (151) to solve the convolution integral in Fourier space, and 60 anti-transformed to obtain a new guess for the strain-rate eld, etc. Further- o o more, assuming (x) = (x) as the initial guess for an auxiliary stress ij ij eld associated with the compatibility constraint, the Augmented Lagrangian procedure (eqs. (142)-(143)) can be readily applied to iteratively determine the local elds. Microstructural changes, including porosity evolution, can then be tracked using the updating algorithms described in Lebensohn et al. (2013). 6. Multiscale modelling of polycrystals Simulation of polycrystals by means of either mean- eld homogenization schemes or full- eld homogenization of an RVE of the microstructure can only consider homogeneous stress or strain boundary conditions. Thus, it cannot be used to determine the mechanical response of polycrystals with arbitrary geometry or subjected to complex mechanical loads. Multiscale frameworks in which the single crystal (microscopic) behavior is connected with the specimen (macroscopic) response through an intermediate (meso- scopic) scale are the most direct solution to link the non-homogeneous defor- mation of an structure with its microstructure. The behavior at mesoscopic scale is obtained from the homogenized response of a polycrystal, that repre- sents the behavior of each (polycrystalline) material point in the macroscale. These models allow considering microstructural e ects on the response of a polycrystalline specimen and they are only applicable to problems with a clear separation of scales, in which the length-scale associated with the gra- dients of the mechanical elds at the macroscale has to be large compared with the length-scale of the polycrystalline microstructure, i.e. the grain or sub-grain size. The models and numerical methods at both ends of this multiscale prob- lem are well-established and readily available. The mechanical behavior of single crystals is modelled using any of the CP models presented in section 3, while the macroscopic behavior at the macroscopic scale, possibly of com- plicated geometry and undergoing complex boundary conditions, is modeled with the FEM using large-strain kinematics. Di erent avours of the multi- scale simulation strategy are determined by the model used at the mesoscopic scale. For example, they can be based on full- eld homogenization, in which the mechanical response at each integration point is given by numerical sim- ulation of a polycrystalline RVE associated to that integration point. Within this group, nite elements have been used to perform this mesoscale compu- 61 tational homogenization for di erent polycrystalline microstructures (Miehe et al., 1999; Werwer and Cornec, 2000; Kouznetsova and Geers, 2008). In the general case of heterogeneous materials (of which a polycrystal is a particu- lar case), these approaches are sometimes referred to as FE models (Feyel, 2003). However, FEM-based homogenization at the mesoscale requires a huge computational e ort for the whole multiscale calculation. An obvious direc- tion for reducing this cost is the use of more ecient full- eld models based on FFT to obtain the average behavior of the RVE in each integration point (section 5.2). A more complex alternative is the use of model order reduction for the mecanical behavior at the microscale. Model order reduction has been mainly used for linear heterogeneous materials or more simple constitutive models at the microscale, but it has been recently applied to polycrystals (Michel and Suquet, 2016). Finally, the computational cost of concurrent multiscale models is strongly reduced if, instead of a numerical model, a mean- eld homogenization approach (section 4) is used at the microscale. In this latter case, the polycrystal response at each integration point is ob- tained from an estimate of its e ective behavior, based a non-space resolved statistical representation of the microstructure and homogenization theory to connect the single crystal and polycrystal levels. Another important aspect of the implementation at mesoscopic level within the multiscale analysis, as discussed by Van Houtte et al. (2006), concerns whether the polycrystal model is embedded in the FE computation (i.e. interrogated on the y each time the behavior of the corresponding material point is needed by the FE analysis) or used in a hierarchical fash- ion (i.e. the FE model utilizes a number of parameters which are identi ed in advance (pre-computed), using the mesoscale model) (Van Houtte et al., 2006, 2009). Embedded models are obviously more accurate to follow the microstructural evolution and how the later a ects the macroscopic behav- ior, while hierarchical approximations are more ecient. In connection with these two di erent strategies, it is also worth mentioning a promising inter- mediate approach, known as the adaptative sampling method (Barton et al., 2008), which consists in embedding a lower length-scale model (e.g. at poly- crystal level) and storing the response obtained for a given microstructure each time this model is interrogated. In this way, a microstructure-response database is populated, which can be used instead of a new call to the poly- crystal model when the current microstructure is not far from an already computed point of this microstructure-response space. 62 Currently, one of the most ecient multiscale approaches is based in the combination of FEM at the macroscopic level with a VPSC approximation at the mesoscale. Along these lines, Segurado et al. (2012) developed a semi- implicit elasto-visco-plastic multiscale framework to simulate the plastic de- formation of polycrystalline structures. This modelling framework have been later extended for di erent type crystals and applied to simulate thermome- chanical processes on full components (Knezevic et al., 2016a; Zecevic et al., 2017a; Knezevic and Beyerlein, 2018). Due to the relevance of this work, the main equations of the framework developed in Segurado et al. (2012) will be reviewed here, including some implementation aspects in a standard nite element code. 6.1. Kinematics The mesoscopic material model based on VPSC (Segurado et al., 2012) was developed using both small-strain and nite-deformation frameworks. For nite deformations, a hypoelastic constitutive formulation was used. Hypoelastic models are a common approach to extend small-strain plas- ticity models to the nite strain range by recasting the original evolution equations in terms of objective stress rates that preserve frame invariance (de Souza Neto et al., 2008). This approach to de ne user materials is a common option in most of commercial codes as ABAQUS, i.e. the one cho- sen for the numerical implementation of the model, due to its simplicity and usability. The application of the hypoelastic framework in this mul- tiscale model allowed using small-strain kinematics inside the VPSC-based user-de ned material subroutine in combination with a co-rotational refer- ence frame that is determined and stored to account for large strains and rotations. In small strains, the total strain increment () is decomposed in elastic ( ) and viscoplastic ( ) parts, such that the constitutive el vp relation at each polycrystalline integration point is given by: = + = L + (155) el vp vp where L is the elastic sti ness tensor of the polycrystalline material point, is the stress increment and = () is computed with the VPSC vp vp model for each polycrystalline material point. For large deformations, the stress increment is based on the interval integration of the Jaumann rate of the Cauchy stress. in consistency with ABAQUS's treatment of rotations. A local coordinate system, de ned by the unit vectors e ( = 1; 3) in which the polycrystal behavior is computed, undergoes rotations that need to 63 be accounted for by the FEM. This system will be referred to as co-rotational system. In the adopted hypoelastic solution of the nite-deformation prob- lem, the material model is formulated incrementally, using the Jaumann rate of the Cauchy stress as the objective rate to de ne the stress increments. In this case, the evolution of the unit vectors of the local coordinate system rotating with the polycrystal is given by: e _ = ! _ e (156) where ! _ is the total rotation-rate of the polycrystalline material point. The integration in a time increment (from time t to t + t) gives an incremental rotation matrix R that rotates the local reference system from t to t + t, t+t t R = RR (157) t+t where R is the rotation matrix from the local coordinate system at t = 0 to its orientation at time t + t, e.g. t+t t+t 0 e = R e (158) Similarly to eq. (155), the velocity-gradient L for large strains is decom- posed into elastic and viscoplastic components L = L + L (159) el vp where: L = _ + ! _ (160) el el and L = _ (161) vp vp Elastic and viscoplastic constitutive laws involve only the symmetric com- ponent of the deformation, while the antisymmetric component is ambiguous. Consequently, the full rotation-rate in eq. (160) was assigned to the elastic part of the velocity-gradient for convenience. The Jaumann rate of the Kircho stress, , which is the stress measure work-conjugate with the rate of deformation, is de ned by = _ ! _ + ! _ : (162) 64 If plastic strains are much larger than elastic ones, deformation can be considered isochoric (J = det (F) = 1 and = J ). Thus, is related to the elastic strain-rate by = L ( _ _ ) : (163) vp The FE solution requires that the constitutive relation in eq. (163) is expressed and integrated in a xed coordinate system. If the coordinate system associated to every polycrystalline material point is unique and given by e , this reference frame is used as the required xed coordinate system. Integrating eq. (163) from t to t + t gives t = L ( ) (164) vp In this formulation, the polycrystal elastic sti ness and the viscoplastic strain increment are calculated by means of the VPSC model in the co- t+t rotational system. Therefore, use is made of R to express these magni- tudes in the xed coordinate system, e.g. if is the viscoplastic strain vp increment computed by means of the VPSC model (the symbol \*" denotes magnitudes expressed in the co-rotational system), in eq. (164) is ob- vp tained as t+t t+t = R R : (165) vp vp In eq. (164), is the logarithmic strain increment, obtained from given by the FE analysis as = log (V) (166) with V being the left Cauchy stretch tensor associated with the deformation gradient increment from t to t + t, such that F = VR. Hence, t+t t+t = R R : (167) Finally, for completeness, once the right-hand of eq. (164) is obtained by computing the left-hand magnitudes in the xed reference system, the stress expressed in co-rotational axes (as needed by VPSC as part of the proposed algorithm, see next section) is obtained as r T t+t t+t t t+t = R + t R : (168) 65 6.2. VPSC-based UMAT implementation The approach at the macroscopic level in the FEM implementation is standard: the applied load is divided in increments, and the equilibrium at each increment is obtained by means of the FEM analysis in an iterative fashion, using a global non-linear solver. The load increment is controlled by the time, and once the problem has being solved at time t, the solution for the next time increment requires the polycrystal model to provide a tangent tg sti ness (Jacobian) matrix L = @/@ for each material point, so the FEM can compute an initial guess for the nodal displacements at t + t. The strain increments obtained from that prediction for each material point, FE t , together with the stress and the set of internal state variables q corresponding to the previous increment, are used inside the UMAT to calculate a new guess for the stress and the Jacobian at t + t. When convergence in stress equilibrium is achieved by the global non-linear scheme, the new values (at t + t) of the stresses, the internal variables, and the Jacobian matrix are accepted for every node, and the calculation advances FE to the next increment. For a given , the VPSC-UMAT is based on the minimization of an ad hoc residual, de ned below. t+t t t = + L = + L ( ) (169) el vp where L is the elastic sti ness of the polycrystal. In the present context, the natural choice for L is to use the elastic self-consistent (ELSC) estimate, given by eq. (97), D E (r) (r) L = L B (170) (r) (r) where L and B are the elastic sti ness and stress localization tensors, eq. (82), of SR grain (r). Note that an ELSC calculation for the determination of L is implemented in the VPSC code at the beginning of each deformation increment. In this way, the texture changes are also accounted for in the determination of the sti ness of the polycrystalline material element. Combining eq. (96) with the viscoplastic constitutive relation coming from VPSC, eq. (93), leads to, 1 (px) t t = L + t _ + ; q : (171) vp i The residual X () can be de ned as a (non-linear) function of the t+t t stress increment = , which the UMAT should return to the 66 FEM, as the constitutive response of the polycrystalline material point to FE the trial strain increment FE 1 (px) t t FE X () = = L + t _ + ; (172) vp i FE The condition X () = 0 (i.e. = ), is obtained using a Newton- Raphson (NR) scheme to solve this non-linear algebraic equation. The cor- responding Jacobian, J ; is given by: NR (px) @ _ @X () vp 1 t t = J () = L + t + ; q = NR @ () @ () 1 (px) t t = L + tM + ; q (173) k 1 Hence, given a guess for the stress increment, the new guess is obtained as: k k 1 1 k 1 k 1 = J X : (174) NR Moreover, using the rst equalities of eqs. (172) and (173), FE @X () @ () J = = = (175) NR @ () @ () @ () the Jacobian matrix that the VPSC-based UMAT should pass to the FEM is obtained, @ () 1 tg 1 1 (px) L = = J = C + t M : (176) NR @ () Eq. (176) provides a closed expression for the FEM Jacobian, as a func- tion of the viscoplastic tangent moduli (which is calculated as part of the VPSC algorithm), the elastic sti ness of the aggregate, and the FE time in- crement. The use of this expression greatly reduces the overall computational cost because the polycrystal stress and the elasto-viscoplastic tangent sti - ness tensor are obtained from the same calculation loop. Moreover, the use of eq. (176) for the FEM Jacobian usually provides quadratic convergence to the macroscopic non-linear equations. Finally, it must be remarked that the above procedure is semi-implicit: implicit in the stress value, because internal equilibrium is checked at the end 67 of the increment, but explicit in the internal variables q (such as grain orien- tations, morphology and hardening variables). The reason for the proposed algorithm to remain explicit in those internal variables is related to the dif- culty and computational cost of deriving a residual and a Jacobian matrix that would contain literally thousands of internal variables. In any case, the explicit updating of the internal variables does not appear to change signi - cantly the convergence of the macroscopic model in most cases, as shown by Segurado et al. (2012). 7. Application to virtual testing of polycrystals In order to demonstrate the potential of CHP, applications of the tools and strategies presented above are illustrated in this section, with exam- ples of simulations of the mechanical behavior of polycrystals under di erent conditions. 7.1. Strength The mechanical behavior of materials is usually characterized through a uniaxial stress-strain curves that contain the information about the most important mechanical parameters: elastic modulus, yield strength, strain hardening rate and ductility. Thus, the prediction of uniaxial stress-strain curves is obviously the rst step to validate the capabilities of CHP and there are many examples in the literature for di erent polycrystalline mate- rials (see, among many others, Roters et al. (2010b,a); Zhang et al. (2015); Cruzado et al. (2015); Ghosh et al. (2016); Shahba and Ghosh (2016)). From a Materials Science perspective, the information required to carry out an ac- curate prediction of the stress-strain curve includes a detailed description of the microstructure, and of the hardening parameters that dictate the plastic deformation of the constituent single crystals, as well. While microstrutural information can be easily obtained using modern characterization techniques, the quanti cation of hardening behavior has to be achieved using di erent strategies based on: a) micromechanical tests, b) inverse optimization, or c) multiscale models. While multiscale approaches to predict the CRSSs of the di erent slip systems in engineering materials are not yet mature, one example of each of the other two strategies is detailed below. 68 7.1.1. Strength of IN718 Ni-based superalloy Polycrystalline IN718 is a Ni-Fe based superalloy widely used for struc- tural applications up to 650-700 C because of its good castability and weld- ability, high strength and corrosion resistance. The microstructure of IN718 is made up by a Ni FCC solid solution which contains a dispersion of nm- 0 00 sized (Ni (Al,Ti)) and (Ni Nb) coherent precipitates within the grains 3 3 together with m-sized metal carbides and phase (Ni Nb) particles at grain boundaries (Fig. 5) (Radavich, 1989; Song et al., 2005). The volume frac- 0 00 tions of and phases are in the range 3-5% and 10-20%, respectively, depending on the bulk alloy composition, the heat-treatment and the degree of element segregation (Fayman, 1987). The strength of the alloy in this case (and in the case of many other Ni-based superalloys) is provided by the 0 00 interaction of the dislocations with the ne distribution of and . The precipitate size and spacing is of the order of 10 - 20 nm in wrought IN718 (Fig. 5b), which stands for the critical length scale that controls the me- chanical behavior and, under these circumstances, micromechanical testing techniques that probe the mechanical response of the material in volumes of several m (much larger than the spacing between precipitates) can be used to determine the CRSS. Examples of these strategies in metallic alloys include micropillar compression (Cruzado et al., 2015) and nanoindentation (Hidalgo-Manrique et al., 2015; S anchez-Mart n et al., 2014b; Zambaldi et al., 2015; Vachhani and Kalidindi, 2015; Patel and Kalidindi, 2017). γ'' phase Carbides δ phase γ' phase 50 µm 20 nm Figure 5: Microstructure of IN718 Ni-based superalloy. (a) Polycrystal grain structure showing the distribution of metal carbides and phase at the grain boundaries, (b) Distri- 0 00 bution of and precipitates within the Ni FCC solid solution. Reprinted from Cruzado et al. (2015). The strategy to predict the mechanical behavior of polycrystalline Ni- 69 based superalloys is depicted in Fig. 6 (Cruzado et al., 2015). The mi- crostructural information included the grain size distribution obtained from optical microscopy images, that can be easily translated into the actual 3-D distribution in the case of equiaxed grains (Heilbronner and Bruhn, 1998), while the grain orientation was provided by X-ray di raction measurements. The cubic RVE of the microstructure was generated from the Voronoi tessel- lation of a set of points, which was obtained from a Monte Carlo algorithm, so the grain size distribution in the RVE followed the experimental log-normal grain size distribution, following the strategy presented in section 2.1. The grain orientations within the RVE were assigned to match the measured ODF. Figure 6: Coupled experimental-simulation strategy to predict the mechanical properties of IN718 Ni-based superalloy using computational homogenization of polycrystals. (Cruzado et al., 2015). Deformation of the Ni solid solution occurred by dislocation slip along the 12 f111g<110> slip systems of the FCC lattice. The crystal was assumed to behave as an elasto-viscoplastic solid in which the plastic slip rate in the slip system k, _ , follows a phenomenological power-law dependency according 70 to (eq. (15)) 1=m k k _ = _ sgn( ) (177) where _ stands for the reference strain rate, g is CRSS of the slip system k at the reference strain rate and m the strain rate sensitivity parameter. The evolution of the CRSS of a given slip system, g , is expressed as, k j g _ = h q _ (178) kj where h stands for the self hardening modulus and q are the interaction kj parameters that stand for the in uence of hardening between di erent slip systems. The evolution of the self hardening was described according to the Voce model, h h h 0 s 0 h ( ) = h + h h + exp (179) s 0 s s 0 s 0 where h is the initial hardening modulus, the initial yield shear stress, 0 0 the saturation yields shear stress, h the saturation hardening modulus at s s large strains and stands for the accumulated shear strain in all slip systems, which is given by eq. (18). The parameters of this phenomenological constitutive equation were ob- tained from compression tests in single crystals carried out on micropillars with circular section extracted by means of focus ion beam milling from crystals oriented in di erent orientations. Initially, micropillars of di erent diameters oriented for single slip along the < 235 > or < 245 > directions were tested and the initial CRSS is plotted as a function of the micropillar diameter in Fig. 7a. It was found the the initial CRSS (and the strain hard- ening rate, see Cruzado et al. (2015)) were independent of the micropillar size (and, thus, representative of the bulk properties of the alloy) when the micropillar diameter was > 3 m. Thus, micropillars of 5 m in diameter were tested at di erent strain rates to determine the strain rate sensitivity parameter m (Fig. 7b). In addition, micropillars oriented for double slip (coplanar and non coplanar) as well as for multiple slip were tested. These experimental curves were used to determine the parameters of the constitu- tive model for IN718 (which are shown in Table 1) by comparison with nite 71 _ m h h q 0 0 s 0 s s (MPa) (MPa) (GPa) (GPa) 10 0.017 466 599 6 0.3 1 Table 1: Parameters of the crystal plasticity model for IN718 Ni-based superalloy at ambient temperature obtained from micropillar compression tests. From Cruzado et al. (2015). element simulations of the micropillar compression tests using the constitu- tive equation presented in eqs. (177) { (179). Figure 7: (a) Initial CRSS as a function of the micropillar diameter for micropillars oriented for single slip along the < 235 > or < 245 > directions. One deformed micropillar showing the slip traces is depicted within the plot. (b) E ect of strain rate on the stress vs. plastic strain curve in micropillars of 5 m in diameter oriented for single slip. From Cruzado et al. (2015). The RVE used to determine the mechanical behavior of polycrystalline IN718 by means of computational homogenization is depicted in Fig. 8a. It contained 210 grains and each grain discretized with at least 610 10-node quadratic tetrahedral elements. The grain size distribution and grain orienta- tion followed the experimental data and simulations under periodic boundary conditions were carried out with four di erent realizations of the grain dis- tribution. The di erences among them were below 1.3% and the average true stress vs. logarithmic strain curve under uniaxial compression obtained 72 from computational homogenization is compared with the experimental data in Fig. 8b. The agreement between them was fairly good: the maximum di erence in the compressive ow stress was below 4% and the prediction of the strain hardening rate was very accurate, validating the approach to predict the mechanical behavior of the alloy. (a) (b) Experiment Simulation 0 0.04 0.08 0.12 0.16 0.2 Compressive strain Figure 8: (a) RVE of polycrsytalline IN718 showing the nite element discretization (b) Experimental and computational homogenization results of the true stress vs. logarithmic strain curve in uniaxial compression of IN718 at ambient temperature. From Cruzado et al. (2015). 7.1.2. Strength of AZ31 Mg alloys Computational homogenization of HCP polycrystals to determine the strength presents several additional diculties, as compared with cubic ma- terials. Firstly, HCP have di erent slip systems (basal, prismatic, pyramidal, etc.) with di erent slip resistances that have to be determined independently. Secondly, they often exhibit deformation twinning, which leads to very di er- ent hardening rates. Finally, most HCP metallic alloys present a very strong size e ect during nanomechanical characterization (either using nanoindenta- tion (S anchez-Mart n et al., 2014b; Eswar Prasad et al., 2014; Nayyeri et al., 2017) or micropillar compression (Byer and Ramesh, 2013; Liu et al., 2017)). Thus, the actual values of the CRSS for slip and twinning in the bulk alloy cannot be obtained using the strategy presented in the previous section. Compressive stress (MPa) RVE microstructure: Comparison with grain size, shape, texture experiments Numerical simulation (FEM; FFT) of mechanical tests exp. sim. C-ND T-RD T-ND T-ND/RD45 0 0.02 0.04 0.06 0.08 0.1 0.12 Strain simulation experiments Single crystal plasticity model Error minimization (machine learning) single crystal Error objective properties function Figure 9: Inverse optimization strategy to predict the mechanical properties of AZ31 Mg alloy using computational homogenization of polycrystals. (Herrera-Solaz et al., 2014b). Under these circumstances, the best option to obtain the actual values of the CRSS for slip and twining is to use inverse optimization strategies, as the one depicted in Fig. 9 for AZ31 Mg alloy (Herrera-Solaz et al., 2014b). The rst step of the computational homogenization is the construction of an RVE which includes all the relevant microstructural details (grain size, shape, orientation) and a CP model that takes into account all the actual de- formation mechanisms, namely basal, prismatic and pyramidal slip, together with tensile twinning. The plastic slip rate for a given slip system k follows a power-law according to eq. (15). Similarly, the twinning rate, f , also follows a viscous law given by eq. (55). The activation of each deformation mode depends on the CRSS (g or g ) for each slip or twinning system. The initial values (in absence of previous plastic deformation) of the CRSSs are given by or for the slip system 0 0 k or the twin system , respectively. A phenomenological hardening model was considered for the evolution of the CRSSs, which is able to reproduce Stress (MPa) the di erent stages of single crystal deformation (Kothary and Anand, 1998). The evolution of the CRSSs, g for slip and g for twinning, are expressed by equations (180) and (181), respectively, 3 a a tw j sl k j tw g _ = q h 1 j _ j + q h 1 j _ j (180) sl sl tw sl 0 0 j tw sat sat j=1 tw tw g _ = q h 1 f (181) tw tw tw tw sat where and stand for the resolved shear stress on the corresponding slip (either basal, prismatic or pyramidal) and twinning system, _ is the plastic shear strain rate along system j, f is the rate of the volume fraction transformation in the twin system , and is the characteristic shear of tw the twinning mode ( = 0.129 for extension twinning of Mg alloys, Herrera- tw Solaz et al. (2014b)). It should be noted that extension twinning is a polar mechanism and it will only take place when the applied deformation leads to extension of the c axis of the HCP lattice. The di erent parameters in these equations de ne the contributions aris- ing from self-hardening and latent hardening. The self-hardening of a given slip (k) or twinning ( ) system is de ned by three terms: the saturation stress, , the initial hardening rate h and the hardening exponent a. The sat 0 latent-hardening contribution to slip due to slip in other systems is intro- duced with the coecient q whereas the contribution induced by twinning sl sl is given by q . The model only takes into account the e ect of twinning tw sl on slip and it is assumed that slip does not in uence twinning (q = 0) sl tw following Capolungo et al. (2009) and Zhang and Joshi (2012). One obvious diculty is to nd the actual values of the many parameters of the model. In most cases, they have been tted by comparison with experimental results on polycrystals using a trial and error approach (Miehe et al., 1999, 2002; Segurado et al., 2012; Tom e et al., 2001b, 2002; Pei et al., 2011; Kalidindi et al., 2010; Wu et al., 2010; Tom e et al., 2001a; Segurado and LLorca, 2013) but nding the optimum parameter set is neither easy nor a unique result is guaranteed. In fact, it is not unusual to nd that di erent authors report di erent (or even contradictory) values of the CRSSs of slip and twinning for similar materials. To overcome this limitation, Herrera- Solaz et al. (2014b, 2015) proposed an optimization algorithm which is based in the construction of an objective error function O( ) de ned as 75 p O( ) = jy f (x ; )j = ky f ( )k : (182) i i i=1 where x ; y are a set of p points de ning an experimental data set of the i i polycrystal (for instance, stress-strain curves in di erent orientations, pole gures after deformation, evolution of the fraction of twinned material, etc.) and y = f (x ; ) are the predictions of the corresponding data obtained i i by means of the numerical simulation of an RVE of the microstructure in which the mechanical behavior of each grain within the polycrystal is de- ned by a CP model, which depends on a set of parameters . The mini- mization of the error function can be carried out using di erent algorithms, such as Levenberg-Marquardt (Herrera-Solaz et al., 2014b, 2015) or Nelder- Mead (Zambaldi et al., 2012; Chakraborty and Eisenlohr, 2017). It should be noted that novel machine-learning strategies are very promising to carry out this optimization process (Mueller et al., 2017) and also that the optimiza- tion strategy often requires several iteration loops. Thus, ecient numerical algorithms based on FFT will become critical within this framework. The ability of the inverse optimization strategy to provide accurate values of the CRSS in complex HCP materials, like rolled AZ31Mg alloys processed by hot rolling, was demonstrated by Herrera-Solaz et al. (2014b), Fig. 10a. The average grain size of the rolled AZ31 Mg alloy was 25 m and the pole gure of the as-rolled material is plotted in Fig. 10b. It shows the strong basal texture typical of rolled Mg alloys, with the c axis parallel to the normal direction (ND). Also, the spread prismatic poles appear along the rolling (RD) and transverse (TD) directions. Specimens for tension and compression experiments along di erent orientations indicated in Fig. 10a were machined from the plate and representative stress-strain curves of the mechanical tests can be found in Fig. 11. These curves show the strong plastic anisotropy of Mg alloys, which is triggered by the limited number of slip systems and by the polar nature of extension twinning, which is only activated when deformation leads to an extension of the c axis. As a result, deformation of wrought Mg alloys is markedly dependent on the orientation, and di erent slip systems (and in di erent order) are activated as a function of the loading direction (either tension or compression). For instance, tensile tests along ND show the typical S-shape of the stress-strain curve characteristic of twinning- dominated deformation, while the compression test along ND was dominated by pyramidal slip. 76 (a) (b) Figure 10: (a) Slab of rolled AZ31 Mg alloy showing the RD, ND and TD orientations as well as the orientation of the sample for the mechanical tests. (b) Pole gure of the rolled AZ31 alloy showing the strong basal texture. Reprinted from Herrera-Solaz et al. (2014b). The stress-strain curves corresponding to the tests along ND in tension and compression and in tension along RD were used to determine the pa- rameters of the CP model using the Levenberg-Marquardt method. Full- eld homogenization was carried out on a RVE of the microstructure dis- cretized with cubic nite elements under periodic boundary conditions. The RVE contained 584 grains and was created with the microstructure generator Dream.3D (2018). On average, each grain was discretized with 7 voxels. The grains were equiaxed and the grain size followed a log-normal distribution with an average grain volume equal to the RVE divided by 584. The exper- imental stress-strain curves in the three di erent cases used to calibrate the constitutive equation (tension and compression along ND and tension along RD) were in very good agreement with the simulation predictions. as shown in Fig. 11a. Moreover, this set of parameters was used to predict the stress- strain behaviour of the specimens tested in tension at 45 from the ND/RD orientations (Fig. 11b). The agreement was again good and demonstrates 77 the capabilities of inverse optimization approach to perform virtual test of polycrystals. 300 300 (a) (b) 250 250 200 200 150 150 100 100 exp. sim. Compression ND 50 50 exp. sim. Tension RD Tension ND Tension ND/RD 45 0 0 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.02 0.04 0.06 0.08 0.1 0.12 Strain Strain Figure 11: (a) Experimental (solid lines) and numerical (broken lines with symbol) stress- strain curves resulting from the optimization procedure. (b) Model prediction of the tensile test in the RD-ND plane at 45 from both orientations. Solid lines correspond to experimental results while broken lines with symbols stand for the numerical simulations. Reprinted from Herrera-Solaz et al. (2014b). Further investigations (Herrera-Solaz et al., 2015) were carried out to an- alyze the robustness of the Levenberg-Marquardt optimization strategy to obtain the CRSSs of the di erent slip modes and extension twinning in an AZ31 Mg alloy. It was found that the results were not in uenced by the properties of the single crystals chosen as starting point for the optimization process in so far they were reasonable for the Mg alloy considered. Very similar values of the CRSSs were obtained regardless of whether the initial plastic behavior of the single crystals was isotropic, anisotropic or highly anisotropic. In addition, it was found that there is a minimum critical in- formation that has to be supplied to the optimization strategy in order to achieve meaningful results. In the case of AZ31 Mg alloy, plastic deformation in the basal plane is controlled by basal slip while plastic deformation along the c axis accommodated by either twinning or a combination of prismatic and pyramidal hc+ai slip. Thus, a minimum of two stress-strain curves (one in which deformation along the c is controlled by twinning and another in which is controlled by prismatic and pyramidal hc+ai slip) are necessary as Stress (MPa) Stress (MPa) input for the inverse optimization algorithm. Including a third input stress- strain curve in the optimization process improved the accuracy of the results from the quantitative viewpoint. 7.2. Fatigue life Fatigue crack nucleation in metallic alloys is a process that involves sev- eral length scales, including the features of dislocation slip and disloca- tion/dislocation interactions at the nm scale and the collective motion of dislocations and the formation of particular dislocation patterns at the mi- croscale. They give rise to the localization of plastic deformation in slip bands within the grains that will be the origin of embryonic cracks (Suresh, 2012). This process is very dependent on the local details of the microstruc- ture (grain size, twins, texture) as well as on the presence of defects (pores or intermetallic inclusions). Thus, the role of the microstructure is much more critical in fatigue that in the case of strength. In the latter, the yield strength and the strain hardening are mainly controlled by the average value of the microstructural features (as given by the average grain size, pole gure, etc.) and little di erences are found in polycrystals with the same average features until localization of damage occurs at the last stages of deformation. On the contrary, fatigue crack nucleation occurs by the localization of damage and mainly depends on the extreme values of the statistical distribution of the microstructural features, rather than on the average values. This leads to a very large experimental scatter among "nominally" identical polycrystals (from the viewpoint of average features) and makes very dicult to establish a solid link between microstructure and properties. Obviously, only simulation techniques that take into account the mi- crostructural details are able to provide the link between microstructure and fatigue properties. This is the reason why microstructure-sensitive computa- tional modelling of fatigue(McDowell and Dunne, 2010; Pineau et al., 2016) is indeed among the applications of CHP that has sparked more interest in the last decade because it provides a unique framework to understand the in uence of the microstructure and of the defects on the fatigue behavior of engineering materials. Moreover, metrics can be develop to predict the fatigue life as a function of the loading conditions and of the microstructure of the polycrystal, providing an engineering tool to relate microstructure and fatigue performance. This opens the path to account for design allowables that are controlled by the inherent microstructural variability in the material and to design optimum microstructures for a given loading condition. 79 Figure 12: Coupled experimental-simulation strategy to predict the fatigue life of poly- crystalline metallic alloys using computational homogenization of polycrystals. The modelling framework to predict the fatigue behavior of polycrys- talline metallic alloys by means of CHP is depicted in Fig. 12. The rst step is to obtain the full micromechanical elds during cyclic deformation by means of the numerical simulation of the mechanical behavior of an RVE. The constitutive model for the single crystals should be able to account for the the main features found during cyclic deformation of metallic alloys, namely Bauschinger e ect, mean stress relaxation, ratcheting and cyclic soft- ening/hardening. This is normally achieved by means of phenomenological models (section 3.1), whose parameters have to be calibrated for each al- loy by comparison with the experimental cyclic stress-strain curves. This task can be done by means of trial-and-error approaches or optimization al- gorithms based in the construction of an objective error function (Cruzado et al., 2017). 80 One critical factor to obtain the parameters of the phenomenological CP model as well as to predict the fatigue life is that the mechanical response of the polycrystal changes with the number of cycles until a stable stress-strain hysteresis loop is attained. Stabilization of the mechanical behavior may take from few of cycles up to a few hundreds or thousands of cycles, depending of the alloy and of the loading conditions (Cruzado et al., 2017). In the latter scenario, the computational cost makes impossible the simulation of the mechanical response of the RVE and di erent cyclic jump strategies have been proposed to reach such large number of cycles without the explicit sim- ulation of each of them (Ghosh and Dimiduk, 2011). Among them, Joseph et al. (2010) proposed a novel wavelet transformation-based multi-time scal- ing (WATMUS) algorithm which allows to reduce the number of simulated fatigue cycles by a factor of up to 100. Other simpler approaches can also be found on the literature. For instance, Cruzado et al. (2017) proposed a methodology, following similar approaches in the context of fretting fatigue (McColl et al., 2004), based on the linear extrapolation of the evolution of the internal variables at the constitutive equation level, followed by one or several stabilization cycles. Once the stabilized cyclic stress-strain hysteresis loop has be determined by mean of computational homogenization of the behavior of the RVE, the prediction of the fatigue life is based in some Fatigue Indicator Parameter (FIP) obtained from the evolution of mechanical elds and internal variables at the local level within the RVE. Di erent FIPs have been reported in the literature to describe the main driving force that controls crack formation and they can be linked, using phenomenological relations, to some stage of the fatigue life of the alloy under study (crack nucleation, small crack propagation or full fatigue life). The most common FIPs found in the literature are related to the accu- mulation of some variables during one fatigue cycle, such as the accumulated plastic strain per cycle, P (x) (McDowell et al., 2003; Manonukul and Dunne, 2004; Sweeney et al., 2012, 2014a, 2015) P (x) = L (x)L (x) dt (183) p p cyc or the maximum plastic shear strain (x) (x) = max j _ (x)j dt (184) cyc 81 which are based exclusively on the local plastic strain elds. The Fatemi- Socie FIP, FS, is given by k k (x) (x) FS(x) = max 1 + k (185) k 2 k k where and stand for the maximum plastic shear strain and the max- imum normal stress perpendicular to the slip plane of the k slip system and is the yield stress. This parameter was originally develop to rationalize the fatigue behavior of laboratory specimens (Fatemi and Socie, 1988, 1989) but it was used as the driving force to nucleate a crack in both Low Cycle Fatigue (LCF) and High Cycle Fatigue (HCF) simulations based on CHP (Bennett and McDowell, 2003; Shenoy et al., 2007; Musinski and McDowell, 2012; Przybyla et al., 2013; Castelluccio and McDowell, 2013, 2014). Other popular FIP is the energy dissipated per cycle, W (x), introduced by Charkaluk et al. (2002); Korsunsky et al. (2007) k k W (x) = max (x) _ (x)dt (186) cyc and the energy density stored in each cycle, G(x) (Sweeney et al., 2012, 2014a, 2015; Wan et al., 2016b), which is expressed as (x) : _ (x) G(x) = p dt (187) SSD GND cyc where is a factor that stands for the fraction of the external work in a cycle that has been accumulated though the creation of dislocation structures. Its value is around 0.05 taking into account the typical fraction of the energy dissipated with respect to the external work (Wan et al., 2016b). Many other FIPs can be found in the literature (Anahid et al., 2011; Sangid et al., 2011; Ozturk et al., 2016) but they are not included here for the sake of brevity. There is no clear rule for selecting one of the di erent FIPs presented to predict the fatigue life and, indeed, di erent authors rely on one or another FIP for similar materials. The FIP that better correlates a fatigue process depends on the length scale considered, the characteristics of the alloy and also on the fatigue stage studied (incubation or small crack growth). Persistent slip bands form driven by the accumulation of plastic strain at the single crystal level and under uniaxial loading and both P and are able to reproduce this observation (Manonukul and Dunne, 2004). 82 The role of stress in the development of persistent slip bands and small cracks is much more important at the polycrystal level and under multiaxial loading and this is better accounted for using other FIPs including both plastic slip and stress. A critical assessment of the ability of the di erent FIPs to correlate crack initiation implies comparison with in situ microscopic experiments (Chen et al., 2018) or with lower scale models (Sangid, 2013). The FIPs de ned in eqs. (183)-(187) are local values computed at each point x of the RVE. This local value depends, however, on the details of the discretization and it is not representative of a fatigue damaged region (Sweeney et al., 2013). From a physical viewpoint, the FIP should be av- eraged over a region representative of the crack incubation zone, while the volume averaging is indeed necessary to avoid spurious stress concentrations and minimize mesh size e ects from the computational perspective, as noted in Castelluccio and McDowell (2015). Three di erent averaging volumes are normally used, namely nite element, grain (Shenoy et al., 2007) or slip band (Castelluccio and McDowell, 2015). In the latter, the FIP is averaged on a slip band of constant thickness parallel to each slip plane k at x. After av- eraging, the maximum value of the averaged FIP in the whole RVE (either the maximum averaged FIP in all the nite elements or in all the grains or in all the slip bands) is used to determine the fatigue life as indicated be- low (Shenoy et al., 2007; Castelluccio and McDowell, 2013; Cruzado et al., 2018b). The last ingredient to predict the fatigue life in the microstructure-sensitive computational fatigue models is to establish the link between the FIPs and the di erent stages of the fatigue life. In this respect, micromechanics-based fatigue models can be roughly classi ed in two groups. In the rst one, dif- ferent relationships are established between the FIPs and the fatigue crack nucleation and propagation stages. In the second one, it is assumed that most of the fatigue life is spent in the nucleation of a crack and a direct relationship is established between the FIPs and the fatigue life. The rst approach was introduced by McDowell et al. (2003) and consid- ers at least three stages in the fatigue life: crack incubation, microstructurally small crack growth (MSC) and long crack growth (LC). The fatigue life N is given by N = N + N + N : (188) N MSC LC where N is the number of fatigue cycles to nucleate a crack at a favourable 83 site in the microstructure whose length is of the order of some microstructural length scale such as the grain size. N stands for the number of cycles MSC spent from the nucleation of the crack until it reaches a length of typically 3-10 times the grain size. Crack nucleation and propagation in these two stages are strongly dependent on the microstructure (grain size, texture, size, shape and spatial disitribution of second phases, etc) (Suresh, 2012). Finally, N is the number of the cycles associated with the propagation LC of a crack which is much longer than the longest critical microstructural length scale (typically the grain size in polycrystals). This stage is not very dependent on the microstructure and can be accounted for with the classical crack propagation theories based on fracture mechanics. Thus, CHP is used in this approach to determine N and N . Di erent FIPs and di erent N MSC relationships between each FIP and either N or N can be established N MSC to take into account that the driving forces for both process do not have to be equivalent. The physical process of crack incubation is not yet fully understood. Thus, the determination of the principal driving forces and the develop- ment of robust initiation criteria based on these forces is based in more or less phenomenological approximations. Plastic deformation is accumulated progressively during cyclic deformation because of irreversible slip leading to the formation of low-energy dislocation arrangements to accommodate the irreversible slip. This results in the strain localization in a small region and in the formation of the so-called persistent slip bands, that are the precursors to crack initiation. The rst relation to link a FIP to the cycles for crack nu- cleation based on dislocation mechanics was proposed by Tanaka and Mura (1981), and it has the form TM N (FIP ) = (189) where the FIP in the original work (Tanaka and Mura, 1981) was the ac- cumulated plastic strain per cycle, eq. (183), d stands for the grain size and K is a material constant. Nevertheless, many other expressions and TM models can be found in the literature (Sangid, 2013) and there is no a clear consensus on the form and of variables that should be included. In most of the models in which the FIPs are explicitly related to N (Shenoy, 2006; Shenoy et al., 2007; Przybyla and McDowell, 2010, 2012; Przybyla et al., 2013; Castelluccio and McDowell, 2013, 2014, 2015), this process is linked with the FIPs related to the plastic strain (P and ) following expressions 84 based on the Tanaka and Mura model. Microstructurally small crack growth is treated within the framework of CHP following traditional expressions of macroscopic fatigue crack growth in which the cyclic crack tip displacement range (CTOD) is used as driv- ing force instead of the stress intensity factor range (McDowell et al., 2003). However, CTOD is dicult to obtain in CHP because it would require to de ne a crack and a ne mesh around it (Tucker et al., 2015) but FIPs computed in standard computational homogenization simulations (without the crack) can be used as an e ective surrogate measure of CTOD. There- fore, the microstructural crack growth rate can be obtained as (Shenoy et al., 2007) da = K FS a (190) CTOD y dN where a is the crack length, an averaged value of the critical resolved shear stresses in all the slip systems, g , FS the Fatemi-Socie FIP, eq. (185), and K a constant. Similar expressions but including a threshold have been CTOD used in Przybyla et al. (2013). Eq. (190) can be integrated from an ini- tial crack length of one grain up to a length in which the crack becomes long with respect to the microstructure. This modelling approach has been widely used for modelling HCF in Ni based superalloys,(Shenoy et al., 2007; Musin- ski and McDowell, 2012; Przybyla et al., 2013; Castelluccio and McDowell, 2013, 2014), Al alloys (McDowell et al., 2003) and Ti alloys (Przybyla and McDowell, 2012). The constant that appears in the models for crack nucleation or short crack propagation presented above has always to be obtained by compari- son with some experimental data. This calibration of the model is, however, problematic because it is very dicult to determine experimentally the actual number of fatigue cycles to nucleate a crack or to propagate the crack in the small crack regime. While this may be feasible by means on in situ fatigue tests coupled with di raction and phase-contrast X-ray microtomography (Birosca et al., 2009; Herbig et al., 2011), they are extremely expensive and time-consuming and the available experimental data are limited to macro- scopic fatigue tests in which the crack nucleation event cannot be detected. A simple approach to overcome this problem was introduced by Manonukul and Dunne (2004) assuming that the fatigue life is mainly dominated by crack nucleation. Thus, the number of cycles necessary for crack nucleation is iden- ti ed with the total fatigue life and the relationship is established between 85 the FIP and the fatigue life. Manonukul and Dunne (2004) also showed that this relationship was characteristic of the polycrystalline C263 Ni-based superalloy and independent of the loading conditions, etc. This direct rela- tionship between the FIP and the fatigue life has been profusely used since then using di erent FIPs and also considering di erent expressions to relate the FIP with the fatigue life. For instance, FIPs based on accumulated plas- tic strain have been used to assess the e ect of elastic anisotropy and other factors in the fatigue life of a ferritic steel (Sweeney et al., 2013). The cyclic strain energy dissipation FIP, W (eq. (186)) was used to study the behaviour under LCF in Inconel 718 superalloy (Cruzado et al., 2018b,a), showing that a microstructure-based fatigue model was able to predict the bilnear Con- Manson regime observed experimentally in this alloy. The same FIP was used to predict the fatigue life of stainless steel medical devices (Sweeney et al., 2012) and of a CoCr alloy (Sweeney et al., 2014a). The stored energy FIP G (eq. (187)) was used within a physically-based SGCP framework to predict the fatigue life in a ferritic steel in (Wan et al., 2016b) and the same modelling approach was also employed in Sweeney et al. (2014b) to study fatigue in a CoCr alloy. The concept of stored energy as driving force for crack initiation has also been proposed by Sangid et al. (2011) for a Ni-based superalloy where the energy stored was computed using a model that includes continuum and atomistic e ects. As a nal remark, it is important to note the statical nature of fatigue predictions based in CHP. The number of grains in the RVE is always much smaller than that in the actual microstructure and therefore the values of the FIPs obtained in the simulation of di erent RVEs corresponding to the same microstructure can be quite di erent. Thus, in general, the results obtained with a single RVE are not fully representative of the fatigue life for a given microstructure and do not provide reliable information about the scatter resulting from the particular orientation arrangement. This limitation is often overcome with the concept of SRVEs (Jeulin et al., 2004) presented in section 2. Fatigue estimations are usually obtained aggregating the results of many di erent RVEs which are statistically representative of the same microstructure (Shenoy, 2006; Castelluccio and McDowell, 2013). Prediction of the fatigue life in IN718 Ni-based superalloy To illustrate the fatigue life prediction process using CHP, an example of prediction of the fatigue life in IN718 Ni-based superalloy in the LCF regime at 400 C is presented next (Cruzado et al., 2017, 2018b,a). The CP model 86 used to represent the cyclic plastic behavior (Cruzado et al., 2017) included the mechanisms leading to Bauschinger e ect, mean-stress relaxation and cyclic softening and a summary of the constitutive equations was already presented in eqs. (25)-(27). The fatigue life estimation was based in the plastic energy dissipated by cycle FIP, eq. (186), averaged in slip bands (Castelluccio and McDowell, 2015). The resulting FIP of a particular RVE, was given by W = max max W (x)dV (191) cyc i cyc i=1;nb k i V i V where k (= 1, 2, 3) corresponds to the three di erent slips systems contained in the slip plane parallel to the band i, V is the volume of that band and nb is the total number of bands in the microstructure, which is 4 times the number of elements in the RVE. The fatigue life estimation for this partic- ular RVE (the number of cycles N ) is obtained as function of W using i cyc a power-law relation, introduced by Cruzado et al. (2018b), to account for the change in deformation mechanism that controls the nucleation of fatigue cracks in this material: from localized deformation in a few grains at small cyclic strain ranges to homogeneous plastic deformation at large cyclic strain ranges. Thus, crit N = (192) (W ) cyc where the parameters W and m were obtained from two independent fa- crit tigue experiments corresponding to two di erent cyclic strain ranges under the same strain ratio R =-1 (Fig. 13a). To this end, the FIPs (W ) ob- " cyc tained by averaging the results of 20 di erent RVEs for each of the two loading conditions selected were used. Finally, the model was used to estimate the fatigue life under di erent loading conditions. Predictions for each loading condition were obtained by averaging the FIPs computed using di erent RVEs. The ability of the proposed model to predict the fatigue life of IN718 alloy at 400 C is depicted in Figs. 13(a) and (b), in which the experimental and the model results for the fatigue life are plotted as a function of the cyclic strain range in tests carried out with R = -1 and 0, respectively. The model was able to predict accurately the fatigue life under small and large cyclic strain ranges 87 for all the loading cases considered. Moreover, the numerical predictions with di erent RVEs led to estimations of the scatter in fatigue life which followed the same trends that the exerpimental observations: the scatter decreased as the applied cyclic strain range increased. 3.5 3.5 (a) (b) 3 3 Model calibration 2.5 2.5 2 2 1.5 1.5 1 1 Experiments Experiments 0.5 (R =-1) 0.5 (R =0) Simulations Simulations " " 0 0 2 3 4 5 2 3 4 5 10 10 10 10 10 10 10 10 Cycles (N ) Cycles (N ) i i Figure 13: Experimental results and model predictions for the fatigue life of IN718 Ni- based superalloy at 400 C. (a) R = -1 and (b) R = 0. The cyclic strain amplitude, , " " is normalized by , the minimum cyclic strain amplitude used in the tests. Reprinted min from Cruzado et al. (2018b). 7.3. Void growth in polycrystals Ductile failure of polycrystalline metals takes place by the nucleation, growth and coalescence of voids (Thomason, 1990). Void nucleation is nor- mally triggered by either fracture or interface decohesion of second phase particles while void growth up to coalescence { that controls the nal duc- tility of the material { is driven by the plastic ow and mainly depends on the stress state and on microstructural details. While the e ect of stress state on void growth and coalescence has been analyzed in detail within the framework of isotropic plasticity models (Benzerga and Leblond, 2010), the in uence of the microstructure (grain size, orientation) is far less understood. Nevertheless, it has been experimentally established that grain anisotropy and grain boundaries play an important role on void growth, particularly for low symmetry crystals (Lecarme et al., 2014; Nemcko and Wilkinson, 2016; Pushkareva et al., 2016). ""="" min ""="" min CHP has been applied to analyze the e ect of the microstructural factors on void growth using the dilatational viscoplastic FFT-based model described at the end of section 5.2. This approach was used to simulate void growth in damaged polycrystalline materials under high stress triaxiality loading that resembles that of typical Cu spall experiments (Escobedo et al., 2011, 2013). The simulation results in two cases in which the matrix behavior was di erent (FCC polycrystal vs. isotropic matrix) were compared with post-mortem orientation images from one of such experiments, showing grain orientations and voids in the region of incipient damage. This provides information on the e ect of the matrix's polycrystallinity on porosity evolution, and identify microstructural e ects on void growth, such as the in uence of the Taylor factor of the crystalline ligaments that link interacting voids. The mechanical behavior of the single crystals in the polycrystal assumed an FCC material with 12 f111g <110> slip systems, a strain rate sensitivity exponent n =5 and a reference strain rate _ = 1, eq. (15). The initial critical resolved shear stress was 100 MPa with a simple linear hardening h kj = 30 MPa, eq. (17). These values are consistent with the ones measured in shock-compressed Cu with 5% pre-strain (Sencer et al., 2005) . Meanwhile, for consistency, the constitutive parameters that describe the behavior of an analogous isotropic matrix material, eq. (148), were n = 5, _ = 1, = 250 0 0 MPa and H = 75 MPa. The polycrystalline RVE included 200 grains with an initial porosity of 1% distributed in 50 spherical voids randomly seeded at grain boundaries. It represents the spall region of a polycrystalline high-purity Cu target fol- lowing compressive shock hardening and subsequent random nucleation of intergranular voids. The isotropic matrix RVE was built keeping the same void distribution and replacing the grains by a homogenous isotropic matrix with the analogous constitutive behavior given above. The initial porosity represented spallation voids that nucleate when in- teracting rarefaction waves produce the required tensile state for damage ini- tiation. From that point on, it was assumed that the material deforms under an axisymmetric stress state, with far- eld stresses > = (where zz xx yy z is the shock direction, normal to the spall plane), and a constant stress- triaxiality X = = , where = ( + 2 ) =3, and = . m eq m zz xx eq zz xx The local elds and the porosity evolution for both the polycrystal and the isotropic matrix cases are plotted in Figure 14 for X = 2.5 . Figs. 14a and c stand for the relative equivalent strain rate elds calculated at the initial and nal deformation steps, respectively, in the case of the polycrystalline 89 FCC matrix. The corresponding results for the isotropic matrix are plotted in Figs. 14b and d. The black arrows indicate the direction of the largest principal stress, normal to the spall plane. The volumetric deformation is fully accommodated by void growth, resulting in a strong localization of (incompressible) plastic deformation in the material surrounding the cavities. This is revealed by the high values (relative to the macroscopic equivalent strain rate) of the plastic strain rate elds, especially between cavities that are close to each other and whose ligaments are relatively well aligned with the direction of largest principal stress. In order to better appreciate the di erences between both cases, 2-D sec- tions from the 3-D simulations of Fig. 14 are plotted in Figs. 15a to c for the porous material with isotropic matrix and in Figs. 15d to f the case of the FCC polycrystalline matrix at di erent levels of porosity. The experimental results of the typical post-mortem EBSD-generated Taylor factor maps of a damaged polycrystalline Cu, from the spall region of the target are also shown in Figs. 15g and h (Escobedo et al., 2011, 2013). Void growth in the case of the isotropic matrix is only a ected by the distance and relative position of the interacting voids. For example, the two pairs of voids within the yellow circle in Fig. 15c interact strongly and grow in a similar way. On the contrary, the voids in the lower pair (marked 2) do interact and grow profusely while the upper pair (marked 1) stop growing, likely due to the presence of a hard crystalline ligament, in the polycrystalline matrix case (Fig. 15f ). Indeed, this ligament is mostly formed by grain 92 in this par- ticular 2-D section, whose initial Taylor factor with respect to tension along the z direction is 3.31, i.e. a hard crystallographic orientation for deform- ing in tension along that axis. Interestingly, the EBSD Taylor factor maps show similar behaviour. Fig. 15g shows two voids that are ideally located to interact, but they did not coalesce, likely because of the hard ligament linking them (note that the orange and red colours represent Taylor factors > 3.1). On the other hand, Fig. 15h shows three voids (the two on the right have already merged) in similar spatial con guration as those of Fig. 15g, but linked by soft orientations (Taylor factor < 2.5), which have coalesced or are in an advanced state of pre-coalescence. It should be emphasized the qualitative character of the above comparison, which is based on two strong assumptions, namely: 1) that the e ect of the 3-D neighbourhood can be neglected, and 2) that the stress state remained axisymmetric throughout the entire deformation process, to justify the use of Taylor factors along the shock (axial) direction. Figure 14: Contour plots of the equivalent strain rate elds during spalling of a porous metal when the matrix is modelled as an FCC polycrystal with intergranular cavities (left), or a homogenous isotropic matrix (right), for X = 2.5. (a-b) Initial, (c-d) nal con gurations. The black arrow indicates the direction of the largest principal stress, normal to the spall plane. Simulations were carried out using full- eld homogenization with FFT. Reprinted from Lebensohn et al. (2013). Figure 15: 2-D sections corresponding to the 3-D simulations in Fig. 14, at di erent levels of porosity, for porous materials with (a-c) isotropic matrix, and (d-f ) FCC polycrystalline matrix. (g-h) Typical post-mortem EBSD-generated Taylor factor maps of dynamically- damaged polycrystalline Cu, recovered from the spall region. Reprinted from Lebensohn et al. (2013). 92 7.4. Multiscale modelling of rolling Multiscale models with homogenization at the mesoscale (described in section 6) are useful to simulate the mechanical behavior of polycrystals sub- jected to non-homogeneous stress states. An example of application of this strategy is the analysis of the spatial variation of the texture evolution during rolling of an FCC polycrystal (Segurado et al., 2012). These variations, in turn, change the local behavior and, as a consequence, determine the redis- tribution of stresses during deformation. This case constitutes a challenge for the proposed approach, since it includes heterogeneous deformation, large strains and rotations, contact and complex local strain histories. The simulation of rolling at the macroscopic level is carried out using the FEM. A parallelepipedic slab of material is rolled between two cylinders leading to 50% reduction in thickness (Fig. 16). The initial aspect ratio of the slab was 4x1 in the rolling-direction/normal-direction plane (RD-ND plane). Plane-strain was imposed in the transverse direction, TD, and symmetry was used to reduce the size of the FEM model. The RD and TD coincided with the x and z axes of the macrocopic FEM. The slab was discretrized with a structured mesh of linear cubes with full integration (C3D8 in Abaqus) using 26 x 9 elements in the RD-ND plane and 1 element in TD. Load was applied in two steps. The initial step consisted in a small vertical displacement of the cylinders in order to establish contact conditions between rolls and plate at the entry and the cylinders rotate at 2 rad/s in the second step until the plate has been completely rolled. The movement of the cylinders was transferred to the plate by friction, and the imposed angular velocity corresponds to an average compression strain-rate component of along the normal direction (ND, y axis). The mechanical behavior of the material at each Gauss point of the FEM was determined using the VPSC approach. The FCC polycrystalline mate- rial was initially represented by 500 randomly-oriented grains which deform plastically in 12 f111g <110> slip systems. The elastic sti ness tensor of the Al FCC crystal was given by L = 108 GPa, L = 62 GPa and L 1111 1212 4444 = 28 GPa. The VPSC constitutive parameters, adjusted to reproduce the experimental behavior of a 6116 Al alloy deformed in uniaxial compression (Tom e et al., 2002) followed a Voce law, eq. (179), characterized by 116 MPa, = 119 MPa, h = 793 MPa and h = 31 MPa, isotropic hardening s o s (q =1, 8k; j), and a strain rate sensitivity exponent n =10. kj The geometry of the rolled plate in the nal stage is shown in Fig. 16, together with the corresponding contour plots of the Von Mises stress, the (a) (b) (c) Figure 16: Multiscale nite element simulation of rolling of an FCC plate, up to 50% reduction. Contour plots of: a) Von Mises stress (in MPa), b) diagonal strain component in ND, . c) shear strain component along RD normal to ND, . Reprinted from yy xy Segurado et al. (2012). 94 J. Segurado et al. / International Journal of Plasticity 28 (2012) 124–140 135 The movement of the cylinders is transferred to the plate by friction, and the imposed angular velocity corresponds to an !1 average compression strain-rate component of 2.1 s along the normal direction (ND, y-axis). The plate of initial size 4 " 1 in the rolling (RD, x-axis)-ND plane, is meshed with a structured mesh of linear cubes with full integration (C3D8 in ABAQUS) using 26 " 9 elements in that plane and 1 element in TD. The geometry in the final stage is shown in Fig. 4, together with the corresponding distributions of Von Mises stress, the diagonal strain component in ND, e , and the shear strain e . Interest- yy xy ingly, the stress does not vanish in the region already laminated and a heterogeneous distribution of residual stresses is left in the material. Both from Fig. 4c and the final shape of the elements it can be concluded that elements near the surface underwent significant shear deformation along RD normal to ND, while elements towards the center were subjected to smal- ler shear distortion. This shear deformation is a consequence of the friction between the metal surface and the cylinders. The different deformation history of the polycrystalline material points through the section determines different evolu- tion of the texture. In order to study this effect, the texture has been tracked in two material points, located near the center and near the surface of the plate, respectively. Both points were taken far from the leading tip in the rolling direction in order to minimize edge effects and obtain results that resemble those of a continuum rolling process. The h111i pole figures of both polycrystalline material points are shown in Fig. 5. It is observed that the textures near the center resemble typical plane-strain compression textures of high-SFE FCC polycrystals, while the texture associated with a point near the surface shows significant rotation with respect to TD of the rolling texture components, consistent with a plane-strain + shear strain history applied to the material point. Similar results were reported by Engler et al. (2000), in terms of experimental Al rolling textures and SA-VPSC texture predictions using FE-generated deformation histories. 4.2. Dimensional changes in bent HCP bars Four-point bending experiments performed on textured Zr (Tomé et al., 2001; Kaschner et al., 2001) and Ti (Nixon et al., 2010a) bars have been used to illustrate the strong effect of texture-induced anisotropy on the mechanical response of HCP polycrystals. These experiments consisted in the bending of bars of square cross-section, cut from the same plate, with their main basal texture component oriented either perpendicular or parallel to the bending plane. This geometry resulted in dif- ferent dimensional changes of the cross-sections of the bars, which can be explained in terms of texture-related ‘‘hard’’ or ‘‘soft’’ responses normal to the fibers along the longitudinal direction. These fibers undergo tension or compression, depend- ing on its vertical position. center surface RD RD TD TD (a) (b) Fig. 5. h11 1i pole figures predicted by ABAQUS–VPSC in the case of rolling of an FCC plate after 50% reduction, for two polycrystalline material points located: (a) near the center, and (b) near the surface of the plate. Lines represent multiple of random distributions (mrd) and dots represent regions with Figure 17: f111g pole gures predicted by multiscale modeling strategy for two polycrys- less than 1 mrd. talline material points located: a) near the center, and b) near the surface of the slab. Lines represent multiple of random distributions (mrd) and dots represent regions with Table 1 VPSC less parameter thans used 1 mrd. in the VPSC–ABAQUS Reprinted simulations fromofSegurado 4-point bendingeof t textu al.red (2012). Zr bars (from Tomé and Lebensohn, 2011). ⁄,PR ⁄,PY ⁄,TT ⁄,CT s (MPa) s (MPa) h (MPa) h (MPa) h h h h n 0 1 0 1 Prhai 14 12 900 40 1 1 10 2 10 Pyrhc + ai 145 192 1684 5 1 1 2 2 10 diagonal strain component in ND, , and the shear strain, . Interestingly, yy xy TT1 90 17 100 30 1 1 10 16 10 CT1 270 30 1000 178 1 1 10 5 10 the stress does not vanish in the region already rolled and a heterogeneous distribution of residual stresses is left in the material. It can be concluded from both Fig. 16c and the nal shape of the elements that elements near the surface underwent signi cant shear deformation along RD normal to ND, while elements towards the center were subjected to smaller shear distortion. This shear deformation is a consequence of the friction between the metal surface and the cylinders. The di erent deformation histories of the polycrystalline material points through the section determine di erent evolution of the texture. This can be appreciated in Figs. 17a and b, which show the <111> pole gures of two material points located near the center and near the surface of the plate, respectively. Both points were taken far from the leading tip in the rolling direction in order to minimize edge e ects and obtain results that resemble those of a continuum rolling process. The textures near the center resemble typical plane-strain compression textures of high stacking fault energy FCC polycrystals, while the texture associated with a point near the surface shows signi cant rotation with respect to TD of the rolling texture components, consistent with a plane-strain + shear strain history applied to the material point. 95 8. Concluding Remarks and Future Developments CHP by means of mean- eld approximations or full- eld simulations of a RVE of the microstructure has become a mature tool in the last twenty years, enabling the simulation of the mechanical behavior of polycrystalline aggre- gates taking into account the e ect of microstructure. This "micromechan- ical" approach is signi cantly di erent from the classical phenomenological plasticity models for polycrystals (for example, J plasticity and associated damage and failure models). These models are known to provide an accu- rate response for the homogenized behavior at the macroscopic level and the parameters are calibrated by means of mechanical tests of polycrystalline specimens. Nevertheless, this classical methodology presents several limita- tions: Each set of parameters is only valid for a given set of testing conditions (temperature, strain rate, loading mode) and microstructural features (grain size distribution and morphology, texture). Any change in these factors re- quires the re-calibration of the model. These phenomenological plasticity models are accurate to predict the homogenized properties (average behavior of many crystals), but fail to de- termine the macroscopic properties when they are controlled by microstruc- tural features (fatigue and creep crack initiation, crack propagation, e ect of defects, etc.). Phenomenological plasticity models assume that the material is homo- geneous and do not include information about the microstructural details, such as grain size and orientation, intragranular distribution of metallurgical phases. Hence, they cannot be coupled with models dedicated to predict the formation and/or evolution of microstructure during thermo-mechanical processing. Obviously, CHP has the potential to overcome these limitations by estab- lishing a clear link between microstructure and mechanical properties that is unfeasible using classical plasticity models. Consequently, CHP is nowadays the standard research tool to analyze the mechanical behavior of polycrys- talline aggregates and is starting to be adopted by industry to solve engi- neering problems (Roters et al., 2010a). This progress in industrial applica- tions has been supported by the availability of di erent open-source software packages that can be used to represent the microstructure of polycrystals (Dream.3D, 2018; Neper, 2018) and to simulate the mechanical behavior us- ing mean- eld (Lebensohn and Tom e, 1993; Lebensohn et al., 2007) and/or 96 full- eld methods based on the FEM or the FFT (Damask, 2018; CAPSUL, 2018). The application of "virtual testing" and " virtual design" method- ologies based on CHP is expected to have a signi cant industrial impact, as it has already happened in the eld of composite materials for aerospace (LLorca et al., 2011; Gonz alez et al., 2017). In particular, it will greatly reduce the number of experimental tests necessary to characterize the me- chanical behavior of new materials and will provide guidelines to modify the microstructure of the latter in order to improve their mechanical performance under di erent loading conditions. In addition, the techniques used in the engineering design will not be based in huge and expensive experimental cam- paigns but on micromechanical models that can account for the variability in the microstructure. Nevertheless, as CHP is used to tackle more and more complex situa- tions, several important limitations of the current strategies to simulate the mechanical performance of polycrystals appear, which should be addressed in the future. The rst important problem to be analyzed is the slip trans- fer across grain boundaries during deformation. There is ample experimental evidence showing that grain boundaries may act as strong barriers for disloca- tions, which are stopped at the boundary leading to the formation of pile-ups, or are easily transferred to the neighbour grain leading to the propagation of deformation on a di erent slip system (Bieler et al., 2014; Bayerschen et al., 2016; Malyar et al., 2017; H emery et al., 2018). Similar e ects have been found for twinning (Hong et al., 2016; Kumar et al., 2016). The nature of the grain boundary from the viewpoint of slip transmission is assumed to depend on three angles that de ne the misorientation between the slip vec- tors, the slip plane normals and the slip plane intersections with the grain boundary plane. Di erent criteria have been proposed but the experimental evidence is still far from being conclusive (Bieler et al., 2014; Bayerschen et al., 2016). Nevertheless, the behavior of grain boundaries may be a deter- minant factor in the mechanical behavior of strongly textured polycrystals (where one type of grain boundary may be dominant) and in the nucleation of damage by shear band formation (Wang et al., 2014) or grain boundary cracking (Boehlert et al., 2012; Cepeda-Jim enez et al., 2015). This is partic- ularly important in the case of low symmetry crystals, in which the limited availability of slip systems may led to large stress concentrations at the grain boundaries if slip transfer is hindered. Moreover, cracks are often nucleated at grain boundaries during fatigue and crack propagation across neighbour grains is critical to determine whether a microstructurally small crack will 97 be arrested or propagate until failure (Musinski and McDowell, 2016; Wan et al., 2016a). Grain boundaries are assumed to be transparent for slip transfer in stan- dard (phenomenological) CP models and, thus, the e ect of the grain bound- ary on the slip transmission is not accounted for. Physically-based CP models can introduce the e ect of grain boundaries through di erent parameters that take into account the increase in dislocation density near the grain boundary (Mayeur et al., 2015; Su et al., 2016; Hauoala et al., 2018). Finally, lower- order SGCP account for the e ect grain boundaries through the density of GNDs generated near the boundaries due to the deformation incompatibil- ity between grains with di erent orientations, while higher-order SGCP can control the dislocation ux across the boundary by means of the higher-order boundary conditions. Thus, either physically-based or SGCP models can be used to simulate the e ect of grain boundaries on the plastic deformation of polycrystals and the development and experimental validation of these strategies will be an important area of research in the near future. Another important development in CHP is the introduction of the "phys- ical" twins in the simulations. So far, simulation of twinning during me- chanical deformation has been taken into account through models that treat twinning as a pseudo-slip system (Kalidindi, 1998; Abdolvand and Daymond, 2013a,b; Chang and Kochmann, 2015). Twinning is triggered by a CRSS on the twin plane and leads to a plastic deformation that depends on the volume fraction of twinned material and on the characteristic shear strain of the twin system considered. Moreover, plastic slip within the twinned regions can also be included. While this approach is accurate from the viewpoint of the ef- fective properties of the polycrystal, it neglects the physical formation of a twinned region which is connected to the parent grain through two coherent grain boundaries. The importance of the presence of the twinned region in the deformation mechanisms have become very clear in the analysis of fatigue crack initiation in FCC metals, which show that fatigue cracks tend to nucle- ate in slip bands on f111g planes parallel and, slightly o set from, annealing twins because of the stress concentration associated with the presence of the twin (Heinz and Neumann, 1990; Stinville et al., 2015, 2016). Obviously, this e ect can only be accounted for if the "physical" twin is included within the grain in the numerical simulation. Annealing twins can be introduced into the initial RVE (although they require very ne meshing of the grains, leading to very large models) (Cerrone et al., 2015; Yeratapally et al., 2016) but twins are often nucleated during deformation and should be introduced 98 on the y during the numerical simulation. This is necessary to account for the hardening e ect of twinning on dislocation slip in a physical way and not through phenomenological models that cannot include the many di er- ent types of interactions. Strategies to achieve this goal have been presented recently (Ardeljan et al., 2015; Cheng and Ghosh, 2017; Liu et al., 2018; Cheng et al., 2018) and they should be further developed. Finally, the crystals in polycrystalline materials subjected to large plastic deformation do not only undergo a change in shape. Grain fragmentation as well as recrystallization processes lead to dramatic microstructural changes during deformation (particularly at high temperatures) that are used during thermo-mechanical processes to manufacture materials with new microstruc- tures and properties. Coupling microstructure evolution with mechanical deformation will allow the "virtual processing" and "virtual design" of novel microstructures with optimized properties for speci c applications and will dramatically increase the add value of CHP. Strategies to achieve these goals have been developed by coupling CHP with phase- eld simulations (Abri- vard et al., 2012a,b; Gvenc et al., 2014; Vondrous et al., 2015; Lim et al., 2016; Admal et al., 2018) or cellular automata (Chuan et al., 2013; Li et al., 2016). They should be developed further to overcome the computational problems associated with the di erent time scales for mechanical deforma- tion and microstructural evolution as well as the experimental and theoretical determination of the driving forces that control the microstructural evolution (grain boundary mobility, nucleation, etc.) under these complex conditions. Acknowledgements This investigation was supported by the European Research Council un- der the European Unions Horizon 2020 research and innovation programme (Advanced Grant VIRMETAL, grant agreement No. 669141), by the Span- ish Ministry of Economy and Competitiveness through the project DPI2015- 67667 and by Los Alamos National Laboratory's Laboratory-Directed Re- search and Development (LDRD) program. References Abdolvand, H., Daymond, M. R., 2013a. Multi-scale modeling and experi- mental study of twin inception and propagation in hexagonal close-packed materials using a crystal plasticity nite element approach. Part I: Average behavior. Journal of the Mechanics and Physics of Solids 61, 783 { 802. 99 Abdolvand, H., Daymond, M. R., 2013b. Multi-scale modeling and experi- mental study of twin inception and propagation in hexagonal close-packed materials using a crystal plasticity nite element approach. Part II: Local behavior. Journal of the Mechanics and Physics of Solids 61, 803 { 818. Abrivard, G., Busso, E., Forest, S., Appolaire, B., 2012a. Phase eld mod- elling of grain boundary motion driven by curvature and stored energy gradients. Part I: Theory and numerical implementation. Philosophical Magazine 92, 3618{3642. Abrivard, G., Busso, E., Forest, S., Appolaire, B., 2012b. Phase eld mod- elling of grain boundary motion driven by curvature and stored energy gradients. Part II: Application to recrystallisation. Philosophical Maga- zine 92, 3643{3664. Acharya, A., Bassani, J., 1995. Incompatible lattice deformations and crystal plasticity. Plastic and Fracture Instabilities in Materials, ASME proceed- ings 57, 75{80. Acharya, A., Bassani, J. L., Beaudin, A., 2003. Geometrically necessary dis- locations, hardening, and a simple gradient theory of crystal plasticity. Scripta Materialia 48, 167{172. Admal, N. C., Po, G., Marian, J., 2018. A uni ed framework for polycrystal plasticity with grain boundary evolution. International Journal of Plastic- ity, In press. Aifantis, E. C., 1987. The physics of plastic deformation. International Jour- nal of Plasticity 3, 211 { 247. Alkemper, J., Voorhees, P. W., 2001. Quantitative serial sectioning analysis. Journal of Microscopy 201, 388{394. Allison, J., Cowles, B., DeLoach, J., Pollock, T., Spanos, G., 2013. Integrated Computational Materials Engineering (ICME): Implementing ICME in the Aerospace, Automotive, and Maritime Industries. The Minerals, Metals & Materials Society, Warrendale, PA 15086. Allison, J., Li, M., Wolverton, C., 2006. Virtual aluminum castings: An industrial application of ICME. JOM 58, 28{35. 100 Anahid, M., Samal, M. K., Ghosh, S., 2011. Dwell fatigue crack nucleation model based on crystal plasticity nite element simulations of polycrys- talline titanium alloys. Journal of the Mechanics and Physics of Solids 59, 2157 { 2176. Anand, L., Kothari, M., 1996. A computational procedure for rate- independent crystal plasticity. Journal of the Mechanics and Physics of Solids 44, 525 { 558. Ardeljan, M., McCabe, R. J., Beyerlein, I. J., Knezevic, M., 2015. Explicit incorporation of deformation twins into crystal plasticity nite element models. Computer Methods in Applied Mechanics and Engineering 295, 396 { 413. Arsenlis, A., Parks, D., 1999. Crystallographic aspects of geometrically- necessary and statistically-stored dislocation density. Acta Materialia 47, 1597 { 1611. Asaro, R., Needleman, A., 1985. Texture development and strain hardening in rate dependent polycrystals. Acta Metallurgica 33, 923{953. Ashby, M. F., 1970. The deformation of plastically non-homogeneous mate- rials. Philosophical Magazine 21, 399{424. Barbe, F., Decker, L., Jeulin, D., Cailletaud, G., 2001a. Intergranular and intragranular behavior of polycrystalline aggregates. Part 1: F.E. model. International Journal of Plasticity 17, 513 { 536. Barbe, F., Forest, S., Cailletaud, G., 2001b. Intergranular and intragranu- lar behavior of polycrystalline aggregates. Part 2: Results. International Journal of Plasticity 17 (4), 537 { 563. Bardella, L., 2006. A deformation theory of strain gradient crystal plastic- ity that accounts for geometrically necessary dislocations. Journal of the Mechanics and Physics of Solids 54, 128 { 160. Bardella, L., Segurado, J., Panteghini, A., Llorca, J., 2013. Latent hardening size e ect in small-scale plasticity. Modelling and Simulation in Materials Science and Engineering 21, 055009. 101 Bargie l, M., Mo scinski, J., 1991. C-language program for the irregular close packing of hard spheres. Computer Physics Communications 64, 183{192. Barton, N. R., Knap, J., Arsenlis, A., Becker, R., Hornung, P. D., Je er- son, D. R., 2008. mbedded polycrystal plasticity and adaptive sampling. International Journal of Plasticity 24, 242{266. Bassani, J., Wu, T., 1991. Latent hardening in single crystals ii. analyti- cal characterization and predictions. Proceedings of the Royal Society of London A 435, 2141. Bayerschen, E., McBride, A. T., Reddy, B. D., B ohlke, T., 2016. Review on slip transmission criteria in experiments and crystal plasticity models. Journal of Materials Science 51, 2243{2258. Bayley, C. J., Brekelmans, W. A. M., Geers, M. G. D., 2006. A comparison of dislocation induced back stress formulations in strain gradient crystal plasticity. International Journal of Solids and Structures 43, 7268 { 7286. Becker, R., 1991. Analysis of texture evolution in channel die compressioni. e ects of grain interaction. Acta Metallurgica et Materialia 39, 1211 { Bennett, V. P., McDowell, D. L., 2003. Polycrystal orientation distribution e ects on microslip in high cycle fatigue. International Journal of Fatigue 25, 27 { 39. Benzerga, A. A., Leblond, J.-B., 2010. Ductile fracture by void growth to coalescence. Advances in Applied Mechanics 44, 169 { 305. Berbenni, S., Taupin, V., Djaka, K. S., Fressengeas, C., 2014. A numer- ical spectral approach for solving elasto-static eld dislocation and g- disclination mechanics. International Journal of Solids and Structures 51, 4157 { 4175. Bertin, N., Tom, C. N., Beyerlein, I. J., Barnett, M. R., Capolungo, L., 2014. On the strength of dislocation interactions and their e ect on latent hardening in pure magnesium. International Journal of Plasticity 62, 72 { 102 Bertin, N., Upadhyay, M. V., Pradalier, C., Capolungo, L., 2015. A t-based formulation for ecient mechanical elds computation in isotropic and anisotropic periodic discrete dislocation dynamics. Modelling and Simula- tion in Materials Science and Engineering 23 (6), 065009. Berveiller, M., Zaoui, A., 1979. An extension to the self-consistent scheme to plastically- owing polycristals. Journal of the Mechanics and Physics of Solids 26, 325{344. Beyerlein, I. J., Tom e, C., 2008. A dislocation-based constitutive law for pure Zr including temperature e ects. International Journal of Plasticity 24, 867{895. Bhattacharyya, A., El-Danaf, E., Kalidindi, S. R., Doherty, R. D., 2001. Evo- lution of grain-scale microstructure during large strain simple compression of polycrystalline aluminum with quasi-columnar grains: Oim measure- ments and numerical simulations. International Journal of Plasticity 17, 861 { 883. Bieler, T. R., Eisenlohr, P., Zhang, C., Phukan, H. J., Crimp, M. A., 2014. Grain boundaries and interfaces in slip transfer. Current Opinion in Solid State and Materials Science 18, 212 { 226. Bilby, B., Crocker, A., 1965. The theory of the crystallography of deformation twinning. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 288, 240{255. Birosca, S., Buere, J.-Y., Garcia-Pastor, F., Karadge, M., Babout, L., Preuss, M., 2009. Three-dimensional characterization of fatigue cracks in ti-6246 using x-ray tomography and electron backscatter di raction. Acta Materialia 57, 5834 { 5847. Bobeth, M., Diener, G., 1987. Static elastic and thermoelastic eld uctua- tions in multiphase composites. Journal of the Mechanics and Physics of Solids 35, 137{149. Boehlert, C., Chen, Z., Gutirrez-Urrutia, I., Llorca, J., Prez-Prado, M., 2012. In situ analysis of the tensile and tensile-creep deformation mechanisms in rolled AZ31. Acta Materialia 60, 1889 { 1904. 103 B ohm, H. J., 2004. A Short Introduction to Continuum Micromechanics. Springer-Verlag Wien. B ohm, H. J., Han, W., 2001. Comparisons between three-dimensional and two-dimensional multi-particle unit cell models for particle reinforced metal matrix composites. Modelling and Simulation in Materials Science and Engineering 9, 47{66. Borg, U., Niordson, C. F., Kysar, J. W., 2008. Size e ects on void growth in single crystals with distributed voids. International Journal of Plasticity 24, 688 { 701. Brenner, R., Beaudoin, A. J., Suquet, P., Acharya, A., 2014. Numerical implementation of static Field Dislocation Mechanics theory for periodic media. Philosophical Magazine 94, 1764 { 1787. Brenner, R., Lebensohn, R. A., Castelnau, O., 2009. Elastic anisotropy and yield surface estimates of polycrystals. International Journal of Solids and Structures 46, 3018{3026. Brisard, S., Dormieux, L., 2010. FFT-based methods for the mechanics of composites: a general variational framework. Computational Materials Sci- ence 49, 663{671. Bronkhorst, C., Kalidindi, S., Anand, L., 1992. Polycrystalline plasticity and the evolution of crystallographic texture in fcc metals. Philosophical Transactions: Physical Sciences and Engineering 341, 443{477. Busso, E. P., Meissonnier, F. T., O'Dowd, N. P., 2000. Gradient-dependent deformation of two-phase single crystals. Journal of the Mechanics and Physics of Solids 48, 2333{2361. Byer, C. M., Ramesh, K. T., 2013. E ects of the initial dislocation density on size e ects in single-crystal magnesium. Acta Materialia 61, 3808 { 3818. Cailletaud, G., 1992. A micromechanical approach to inelastic behaviour of metals. International Journal of Plasticity 8, 55 { 73. Capolungo, L., Beyerlein, I., Tom e, C., 2009. Slip-assisted twin growth in hexagonal close-packed metals. Scripta Materialia 60, 32 { 35. 104 CAPSUL, 2018. https://materials.imdea.org/research/ simulation-tools/capsul/. Castelluccio, G. M., McDowell, D. L., 2013. A mesoscale approach for growth of 3D microstructurally small fatigue cracks in polycrystals. International Journal of Damage Mechanics 23, 791 { 818. Castelluccio, G. M., McDowell, D. L., 2014. Mesoscale modeling of mi- crostructurally small fatigue cracks in metallic polycrystals. Materials Sci- ence and Engineering: A 598, 34 { 55. Castelluccio, G. M., McDowell, D. L., 2015. Microstructure and mesh sen- sitivities of mesoscale surrogate driving force measures for transgranular fatigue cracks in polycrystals. Materials Science and Engineering: A 639, 626 { 639. Cepeda-Jim enez, C. M., Molina-Aldareguia, J. M., P erez-Prado, M. T., 2015. E ect of grain size on slip activity in pure magnesium polycrystals. Acta Materialia 84, 443 { 456. Cermelli, P., Gurtin, M. E., 2001. On the characterization of geometrically necessary dislocations in nite plasticity. Journal of the Mechanics and Physics of Solids 49, 1539 { 1568. Cerrone, A., Stein, C., Pokharel, R., He eran, C., Lind, J., Tucker, H., Suter, R., Rollett, A., Ingra ea, A., 2015. Implementation and veri cation of a microstructure-based capability for modeling microcrack nucleation in lshr at room temperature. Modelling and Simulation in Materials Science and Engineering 23, 035006. Chakraborty, A., Eisenlohr, P., 2017. Evaluation of an inverse methodol- ogy for estimating constitutive parameters in face-centered cubic materials from single crystal indentations. European Journal of Mechanics - A/Solids 66, 114 { 124. Chang, Y., Kochmann, D. M., 2015. A variational constitutive model for slip-twinning interactions in hcp metals: Application to single- and poly- crystalline magnesium. International Journal of Plasticity 73, 39 { 61. 105 Charkaluk, E., Bignonnet, A., Constantinescu, A., Dang Van, K., 2002. Fa- tigue design of structures under thermomechanical loadings. Fatigue and Fracture of Engineering Materials and Structures 25, 1199{1206. Chen, B., Jiang, J., Dunne, F. P., 2018. Is stored energy density the primary meso-scale mechanistic driver for fatigue crack nucleation? International Journal of Plasticity 101, 213 { 229. Cheng, J., Ghosh, S., 2017. Crystal plasticity nite element modeling of discrete twin evolution in polycrystalline magnesium. Journal of the Me- chanics and Physics of Solids 99, 512 { 538. Cheng, J., Shen, J., Mishra, R. K., Ghosh, S., 2018. Discrete twin evolution in mg alloys using a novel crystal plasticity nite element model. Acta Materialia 149, 142 { 153. Cheong, K.-S., Busso, E. P., 2004. Discrete dislocation density modelling of single phase fcc polycrystal aggregates. Acta Materialia 52, 5665{5675. Cheong, K. S., Busso, E. P., Arsenlis, A., 2005. A study of microstructural length scale e ects on the behaviour of fcc polycrystals using strain gradi- ent concepts. International Journal of Plasticity 21, 1797 { 1814. Christian, J., Mahajan, S., 1995. Deformation twinning. Progress in Materi- als Science 39, 1 { 157. Chuan, W., He, Y., Wei, L. H., 2013. Modeling of discontinuous dynamic recrystallization of a near- titanium alloy IMI834 during isothermal hot compression by combining a cellular automaton model with a crystal plas- ticity nite element method. Computational Materials Science 79, 944 { Cruzado, A., Gan, B., Jimnez, M., Barba, D., Ostolaza, K., Linaza, A., Molina-Aldaregu a, J. M., LLorca, J., Segurado, J., 2015. Multiscale mod- eling of the mechanical behavior of IN718 superalloy based on micropillar compression and computational homogenization. Acta Materialia 98, 242 { 253. Cruzado, A., LLorca, J., Segurado, J., 2017. Modeling cyclic deformation of inconel 718 superalloy by means of crystal plasticity and computational 106 homogenization. International Journal of Solids and Structures 122-123, 148 { 161. Cruzado, A., Lucarini, S., LLorca, J., Segurado, J., 2018a. Crystal-plasticity simulation of the e ect grain size on the fatigue behavior of polycrystalline Inconel 718. International Journal of Fatigue 113, 236{245. Cruzado, A., Lucarini, S., LLorca, J., Segurado, J., 2018b. Microstructure- based fatigue life model of metallic alloys with bilinear con-manson be- havior. International Journal of Fatigue 107, 40{48. Cuitino, ~ A., Ortiz, M., 1992. Material-independent method for extending stress update algorithms from small-strain plasticity to nite plasticity with multiplicative kinematics. Engineering Computations 9, 437{451. Dai, H., Parks, D. M., 1997. Geometrically-necessary dislocation density and scale-dependent crystal plasticity. In: Khan, A. (Ed.), Proceedings of Sixth International Symposium on Plasticity. Gordon & Breach, pp. 17{18. Damask, 2018. https://damask.mpie.de. de Botton, G., Ponte-Castaneda, ~ P., 1995. Variational estimates for the creep behavior of polycrystals. Proceedings of the Royal Society of London A 448, 121{142. de Geus, T., Vondrejc, J., Zeman, J., Peerlings, R., Geers, M., 2017. Fi- nite strain t-based non-linear solvers made simple. Computer Methods in Applied Mechanics and Engineering 318, 412 { 430. de Sansal, C., Devincre, B., Kubin, L., 2010. Grain size strengthening in microcrystalline copper: A three-dimensional dislocation dynamics simu- lation. Key Engineering Materials 423, 25{32. de Souza Neto, E., Peric, D., Owen, D. R. J., 2008. Computational Methods for Plasticity. John Wiley and Sons, Ltd. Delaire, F., Raphanel, J. L., Rey, C., 2000. Plastic heterogeneities of a copper multicrystal deformed in uniaxial tension: experimental study and nite element simulations. Acta Materialia 48, 1075 { 1087. Devincre, B., Hoc, T., Kubin, L., 2008. Dislocation mean free paths and strain hardening of crystals. Science 320, 1745{1748. 107 Do sk a r, M., Nov ak, J., Zeman, J., 2014. A periodic compression and re- construction of real-world material systems based on wang tiles. Physical Review E 90, 062118. Dream.3D, 2018. http://dream3d.bluequartz.net. Dunne, F., Rugg, D., Walker, A., 2007. Length scale-dependent, elastically anisotropic, physically-based hcp crystal plasticity: Application to cold- dwell fatigue in Ti alloys. International Journal of Plasticity 23, 1061 { Echlin, M. P., Husseini, N. S., Nees, J. A., Pollock, T. M., 2011. A new femtosecond laser-based tomography technique for multiphase materials. Advanced Materials 23, 2339 { 2342. Eisenlohr, P., Diehl, M., Lebensohn, R. A., Roters, F., 2013. A spectral method solution to crystal elasto-viscoplasticity at nite strains. Interna- tional Journal of Plasticity 46, 37{53. El-Awady, J. A., 2015. Unraveling the physics of size-dependent dislocation- mediated plasticity. Nature Communications 6, 5926. Escobedo, J. P., Cerreta, E. K., Dennis-Koller, D., Patterson, B. M., Bronkhorst, C., Hansen, B., Tonks, D., Lebensohn, R. A., 2011. E ects of grain size and boundary structure on the dynamic tensile response of copper. Journal of Applied Physics 110, 033513. Escobedo, J. P., Cerreta, E. K., Dennis-Koller, D., Trujillo, C. P., Bronkhorst, C., 2013. In uence of boundary structure and near neigh- bor crystallographic orientation on the dynamic damage evolution during shock loading. Philosophical Magazine 93, 833 { 846. Eshelby, J., 1957. The determination of the elastic eld of an ellipsoidal inclusion and related problems. Proceedings of the Royal Society of London A 252, 561{569. Estrin, Y., Mecking, H., 1984. A uni ed phenomenological description of work hardening and creep based on one-parameter models. Acta Metallur- gica 32, 57 { 70. 108 Eswar Prasad, K., Rajesh, K., Ramamurty, U., 2014. Micropillar and macropillar compression responses of magnesium single crystals oriented for single slip or extension twinning. Acta Materialia 65, 316 { 325. Evers, L. P., Brekelmans, W. A. M., Geers, M. G. D., 2004. Non-local crystal plasticity model with intrinsic ssd and gnd e ects. Journal of the Mechanics and Physics of Solids 52, 2379 { 2401. Eyer, D. J., Milton, G. W., 1999. A fast numerical scheme for computing the response of composites using grid re nement. The European Physical Journal Applied Physics 6, 41{47. Fatemi, A., Socie, D. F., 1988. A critical plane approach to multiaxial fatigue damage including out-of-phase loading. Fatigue and Fracture of Engineer- ing Materials and Structures 11, 149{165. Fatemi, A., Socie, D. F., 1989. Multiaxial Fatigue: Damage Mechanisms and Life Predictions. Vol. 159. Springer. Fayman, Y. C., 1987. Microstructural characterization and element partition- ing in a direct aged superalloy (DA718). Materials Science and Engineering A 92, 159{171. Fern andez, A., J erusalem, A., Guti errez-Urrutia, I., P erez-Prado, M., 2013. Three-dimensional investigation of the grain boundary-twin interactions in a Mg AZ31 alloy by electron backscatter di raction and continuum modeling. Acta Materialia 61, 7679{7692. Feyel, F., 2003. A multilevel nite element method (FE2) to describe the response of highly non-linear structures using generalized continua. Com- puational Methods in Applied Mechanics and Engineering 192, 3233{3244. Fleck, N., MUller, G., Ashby, M., Hutchinson, J., 1994. Strain gradient plas- ticity: theory and experiment. Acta Metallurgica et Materialia 42, 475{ Franciosi, P., Zaoui, A., 1982. Multislip in f.c.c. crystals a theoretical ap- proach compared with experimental data. Acta Metallurgica 30, 1627 { 109 Fullwood, D. T., Niezgoda, S. R., Kalidindi, S. R., 2008. Microstructure re- constructions from 2-point statistics using phase-recovery algorithms. Acta Materialia 56, 942 { 948. Geers, M., Cottura, M., Appolaire, B., Busso, E. P., Forest, S., Villani, A., 2014. Coupled glide-climb di usion-enhanced crystal plasticity. Journal of the Mechanics and Physics of Solids 70, 136 { 153. Ghosh, S., Dimiduk, D., 2011. Computational Methods for Microstructure Property Relationship. Springer. Ghosh, S., Shahba, A., Tu, X., Huskins, E. L., Schuster, B. E., 2016. Crystal plasticity FE modeling of Ti alloys for a range of strain-rates. Part II: Image-based model with experimental validation. International Journal of Plasticity 87, 69 { 85. Gong, J., Wilkinson, A. J., 2011. A microcantilever investigation of size e ect, solid-solution strengthening and second-phase strengthening for a prism slip in alpha-ti. Acta Materialia 59, 5970 { 5981. Gonz alez, C., LLorca, J., 2007. Mechanical behavior of unidirectional ber- reinforced polymers under transverse compression: microscopic mecha- nisms and modeling. Composites Science and Technology 67, 2795{2806. Gonz alez, C., Segurado, J., LLorca, J., 2004. Numerical simulation of elasto- plastic deformation of composites: Evolution of stress micro elds and implications for homogenization models. Journal of the Mechanics and Physics of Solids 52, 1573{1593. Gonz alez, C., Vilatela, J. J., Molina-Aldaregu a, J. M., Lopes, C. S., LLorca, J., 2017. Structural composites for multifunctional applications: current challenges and future trends. Progress in Materials Science 89, 194{251. Graham, J. T., Rollett, A. D., Lesar, R., 2016. Fast fourier transform discrete dislocation dynamics. Modelling and Simulation in Materials Science and Engineering 24, 085005. Grennerat, F., Montagnat, M., Castelnau, O., Vacher, P., Moulinec, H., Su- quet, P., Duval, P., 2012. Experimental characterization of the intragranu- lar strain eld in columnar ice during transient creep. Acta Materialia 60, 3655{3666. 110 Gurtin, M. E., 2002. A gradient theory of single-crystal viscoplasticity that accounts for geometrically necessary dislocations. Journal of the Mechanics and Physics of Solids 50, 5 { 32. Gurtin, M. E., 2008. A nite-deformation, gradient theory of single-crystal plasticity with free energy dependent on densities of geometrically neces- sary dislocations. International Journal of Plasticity 24, 702 { 725. Gurtin, M. E., Needleman, A., 2005. Boundary conditions in small- deformation, single-crystal plasticity that account for the burgers vector. Journal of the Mechanics and Physics of Solids 53, 1 { 31. Gvenc, O., Bambach, M., Hirt, G., 2014. Coupling of crystal plasticity nite element and phase eld methods for the prediction of srx kinetics after hot working. Steel Research International 85, 999{1009. Hall, E. O., 1951. The deformation and ageing of mild steel: III Discussion of results. Proceedings of the Physical Society Section B 64, 747{753. Han, C., Ma, A., Roters, F., Raabe, D., 2007. A nite element approach with patch projection for strain gradient plasticity formulations. International Journal of Plasticity 23, 690{710. Han, C.-S., Gao, H., Huang, Y., Nix, W. D., 2005a. Mechanism-based strain gradient crystal plasticity i. theory. Journal of the Mechanics and Physics of Solids 53, 1188 { 1203. Han, C.-S., Gao, H., Huang, Y., Nix, W. D., 2005b. Mechanism-based strain gradient crystal plasticity ii. analysis. Journal of the Mechanics and Physics of Solids 53, 1204 { 1222. Hashin, Z., Shtrikman, S., 1963. A variational approach to the theory of elastic behavior of multiphase materials. Journal of the Mechanics and Physics of Solids 11, 127{140. Hasija, V., Ghosh, S., Mills, M. J., Joseph, D. S., 2003. Deformation and creep modeling in polycrystalline ti6al alloys. Acta Materialia 51, 4533 { Hauoala, S., Segurado, J., LLorca, J., 2018. An analysis of the in uence of grain size on the strength of fcc polycrystals by means of computational homogenization. Acta Materialia 148, 72{85. 111 Hazanov, S., Huet, C., 1994. Order relationships for boundary condition e ects in heterogeneous bodies smaller that the representative volume. Journal of the Mechanics and Physics of Solids 42, 1995{2011. Heilbronner, R., Bruhn, D., 1998. The in uence of three-dimensional grain size distributions on the rheology of polyphase rocks. Journal of Structural Geology 20, 695{707. Heinz, A., Neumann, P., 1990. Crack initiation during high cycle fatigue of an austenitic steel. Acta Metallurgica et Materialia 38, 1933 { 1940. H emery, S., Nizou, P., Villechaise, P., 2018. In situ sem investigation of slip transfer in Ti-6Al-4V: E ect of applied stress. Materials Science and Engineering: A 709, 277 { 284. Herbig, M., King, A., Reischig, P., Proudhon, H., Lauridsen, E. M., Marrow, J., Bure, J.-Y., Ludwig, W., 2011. 3-D growth of a short fatigue crack within a polycrystalline microstructure studied using combined di raction and phase-contrast X-ray tomography. Acta Materialia 59, 590 { 601. Herrera-Solaz, V., Hidalgo-Manrique, P., P erez-Prado, M. T., Letzig, D., LLorca, J., Segurado, J., 2014a. E ect of rare earth additions on the critical resolved shear stresses of magnesium alloys. Materials Letters 128, 199 { Herrera-Solaz, V., LLorca, J., Dogan, E., Karaman, I., Segurado, J., 2014b. An inverse optimization strategy to determine single crystal mechanical behavior from polycrystal tests: Application to AZ31 Mg alloy. Interna- tional Journal of Plasticity 57, 1 { 15. Herrera-Solaz, V., Segurado, J., LLorca, J., 2015. On the robustness of an inverse optimization approach based on the Levenberg-Marquardt method for the mechanical behavior of polycrystals. European Journal of Mechan- ics/A 53, 220 { 228. Hershley, A. V., 1954. The elasticity of an isotropic aggregate of anisotropic cubic crystal. Journal of Applied Mechanics 21, 236{240. Hidalgo-Manrique, P., Herrera-Solaz, V., Segurado, J., LLorca, J., G alvez, F., Ruano, O. A., Yi, S., P erez-Prado, M. T., 2015. Origin of the reversed 112 yield asymmetry in Mg-rare earth alloys at high temperature. Acta Mate- rialia 92, 265{277. Hill, R., 1965a. Continuum micro-mechanics of elastoplastic polycrystals. Journal of the Mechanics and Physics of Solids 13, 89{101. Hill, R., 1965b. A self consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids 13, 213{222. Hill, R., 1966. Generalized constitutive relations for incremental deformation of metal crystals by multislip. Journal of the Mechanics and Physics of Solids 14, 95{102. Hill, R., Rice, J., 1972. Constitutive analysis of elastic-plastic crystals at arbitrary strain. Journal of the Mechanics and Physics of Solids 20, 401{ Hirth, J., Lothe, J., 1982. Theory of Dislocations. Krieger Publishing Com- pany. Hong, X., Godfrey, A., Liu, W., 2016. Challenges in the prediction of twin transmission at grain boundaries in a magnesium alloy. Scripta Materialia 123, 77 { 80. Houtte, P., 1978. Simulation of the rolling and shear texture of brass by the taylor theory adapted for mechanical twinning. Acta Metallurgica 26, 591 { 604. Hu, Z., Rauch, E. F., Teodosiu, C., 1992. Work-hardening behavior of mild steel under stress reversal at large strains. International Journal of Plas- ticity 8, 839 { 856. Hutchinson, J., 1976. Bounds and self-consistent estimates for creep of poly- crystalline materials. Proceedings of the Royal Society of London A 348, 101{127. Idiart, M. I., Moulinec, H., Ponte Castaneda, ~ P., 2006. Macroscopic behavior and eld uctuations in viscoplastic composites: Second-order estimates versus full- eld simulations. Journal of the Mechanics and Physics of Solids 54, 1029{1063. 113 Jeulin, D., Kanit, T., Forest, S., 2004. Representative volume element: A statistical point of view. In: Bergman, D. J., Inan, E. (Eds.), Continuum Models and Discrete Systems. NATO Science Book Series. Springer, pp. 21{27. J.Jung, Yoon, J. I., Kim, J., Latypov, M., Kim, J., H.S, K., 2017. Continuum understanding of twin formation near grain boundaries of fcc metals with low stacking fault energy. npj Computational Materials 3, 21. Joseph, D. S., Chakraborty, P., Ghosh, S., 2010. Wavelet transformation based multi-time scaling method for crystal plasticity FE simulations under cyclic loading. Computer Methods in Applied Mechanics and Engineering 199, 2177 { 2194. Kabel, M., B ohlke, T., Schneider, M., 2014. Ecient xed point and newton|krylov solvers for t-based homogenization of elasticity at large deformations. Computational Mechanics 54, 1497{1514. Kalidindi, S., Bronkhorst, C., Anand, L., 1992. Crystallographic texture evo- lution in bulk deformation processing of fcc metals. Journal of the Mechan- ics and Physics of Solids 40, 537{569. Kalidindi, S. R., 1998. Incorporation of deformation twinning in crystal plas- ticity models. Journal of the Mechanics and Physics of Solids 46, 267{271 and 273{290. Kalidindi, S. R., Knezavic, M., Levinson, A., Harris, R., Mishra, R. J., Do- herty, R. D., 2010. Deformation twinning in AZ31: in uence on strain hardening and texture evolution. Acta Materialia 58, 6230{6242. Khan, R., Zahedi, F. I., Siddiqui, A. K., 2016. Numerical modeling of twin- ning induced plasticity in austenite based advanced high strength steels. Procedia Manufacturing 5, 772{ { 786. Kiener, D., Grosinger, W., Dehm, G., Pippan, R., 2008. A further step to- wards an understanding of size-dependent crystal plasticity: In situ tension experiments of miniaturized single-crystal copper samples. Acta Materialia 56, 580 { 592. 114 Knezevic, M., Beyerlein, I. J., 2018. Multiscale modeling of microstructure- property relationships of polycrystalline metals during thermo-mechanical deformation. Advanced Engineering Materials, 1700956. Knezevic, M., Crapps, J., Beyerlein, I. J., Coughlin, D. R., Clarke, K. D., McCabe, R. J., 2016a. Anisotropic modeling of structural components us- ing embedded crystal plasticity constructive laws within nite elements. International Journal of Mechanical Sciences 105, 227 { 238. Knezevic, M., Zecevic, M., Beyerlein, I. J., Lebensohn, R. A., 2016b. A numerical procedure enabling accurate descriptions of strain rate-sensitive ow of polycrystals within crystal visco-plasticity theory. Computational Methods in Applied Mechanics and Engineering 308, 468{482. Kocks, W., 1975. Thermodynamics and kinetics of slip. Progress in Materials Science 19, 1 { 291. Konrad, J., Zae erer, S., Raabe, D., 2006. Investigation of orientation gradi- ents around a hard laves particle in a warm-rolled fe3al-based alloy using a 3d ebsd- b technique. Acta Materialia 54, 1369 { 1380. Korsunsky, A. M., Dini, D., Dunne, F. P., Walsh, M. J., 2007. Comparative assessment of dissipated energy and other fatigue criteria. International Journal of Fatigue 29, 1990 { 1995. Kothari, M., Anand, L., 1998. Elasto-viscoplastic constitutive equations for polycrystalline metals: Application to tantalum. Journal of the Mechanics and Physics of Solids 46, 51 { 83. Kothary, M., Anand, L., 1998. Elasto-viscoplastic constitutive equations for polycrystalline metals: Application to tantalum. Journal of the Mechanics and Physics of Solids 46, 51{ 67,69 { 83. Kouznetsova, V. G., Geers, M. G. D., 2008. A multi-scale model of marten- sitic transformation plasticity. Mechanics of Materials 40, 641{657. Kowalczyk-Gajewska, K., 2010. Modelling of texture evolution in metals ac- counting for lattice reorientation due to twinning. European Journal of Mechanics - A/Solids 29, 28 { 41. 115 Kreher, W., 1990. Residual-stresses and stored elastic energy of composites and polycrystals. Journal of the Mechanics and Physics of Solids 38, 115{ Kr oner, E., 1958. Berechnung der Elastischen Konstanten des Vielkristalls aus den Konstanten des Einkristalls. Zeitung Physik 51, 504{518. Kumar, M. A., Beyerlein, I. J., McCabe, R. J., Tom e, C. N., 2016. Grain neighbour e ects on twin transmission in hexagonal close packed materials. Nature Communications 7, 13826. Lahellec, N., Michel, J. C., Moulinec, H., Suquet, P., 2001. Analysis of in- homogeneous materials at large strains using fast Fourier transforms. In: Miehe, C. (Ed.), IUTAM Symposium on computational mechanics of solids materials. Kluwer Academic, pp. 247{268. Lautensack, C., 2007. Random Laguerre tessellations. PhD thesis, Universitt Karlsruhe. Laws, N., 1973. On the thermostatics of composite materials. Journal of the Mechanics and Physics of Solids 21, 9{17. Lebensohn, R. A., 2001. N-site modelling of a 3d viscoplastic polycrystal using Fast Fourier Transform. Acta Materialia 49, 2723{2737. Lebensohn, R. A., Brenner, R., Castelnau, O., Rollet, A. D., 2008. Orienta- tion image-based micromechanical modelling of subgrain texture evolution in polycrystalline copper. Acta Materialia 56, 3914{3926. Lebensohn, R. A., Canova, G. R., 1997. Selfconsistent approach for modelling texture development of two-phase polycrystals: application to Titanium alloys. Acta Materialia 45, 3687 { 3694. Lebensohn, R. A., Castelnau, O., Brenner, R., Gilormini, P., 2005. Study of the antiplane deformation of linear 2-D polycrystals with di erent mi- crostructures. International Journal of Solids and Structures 42, 5441{ Lebensohn, R. A., Escobedo, J. P., Cerreta, E. K., Dennis-Koller, D., Bronkhorst, C. A., Bingert, J., 2013. Modelling void growth in polycrys- talline materials. Acta Materialia 61, 6918{6932. 116 Lebensohn, R. A., Hartley, C. S., Tom e, C., Castelnau, O., 2010. Modelling the mechanical response of polycrystals deforming by climb and glide. Philosophical Magazine 90, 567{583. Lebensohn, R. A., Idiart, M. I., Ponte-Castaneda, ~ P., Vincent, P. G., 2011. Dilatational viscoplasticity of polycrystalline solids with intergranular cav- ities. Philosophical Magazine 91, 3038{3067. Lebensohn, R. A., Kanjarla, A. K., Eisenlohr, P., 2012. An elasto-viscoplastic formulation based on fast fourier transforms for the prediction of microme- chanical elds in polycrystalline materials. International Journal of Plas- ticity 32-33, 59{69. Lebensohn, R. A., Montagnat, M., Mansuy, P., Duval, P., Meysonnier, J., Philip, A., 2009. Modeling viscoplastic behavior and heterogenous in- tracrystalline deformation of columnar ice polycrystals. Acta Materialia 57, 1405{1415. Lebensohn, R. A., Needleman, A., 2016. Numerical implementation of non- local polycrystal plasticity using Fast Fourier Transforms. Journal of the Mechanics and Physics of Solids 97, 333{351. Lebensohn, R. A., Tom e, C. N., 1993. A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of poly- crystals: Application to zirconium alloys. Acta Metallurgica et Materialia 41, 2611 { 2624. Lebensohn, R. A., Tom e, C. N., Ponte-Castaneda, ~ P., 2007. Self-consistent modeling of the mechanical behavior of viscoplastic polycrystals incorpo- rating intragranular eld uctuations. Philosophical Magazine 87, 4287 { Lebensohn, R. A., Uhlenhut, H., Hartig, C., Mecking, H., 1998. Mechanical behavior gamma-TiAl-based polysinthetically twinned crystals: microme- chanical modelling and experimental validation. Acta Materialia 46, 4701 { 4709. Lebensohn, R. A., Zecevic, M., Knezevic, M., McCabe, R. J., 2016. Average intragranular misorientation trends in polycrystalline materials predicted by a viscoplastic self-consistent approach. Acta Materialia 104, 228{236. 117 Lecarme, L., Maire, E., Kumar, K. C. A., De Vleeschouwer, C., Jacques, L., Simar, A., Pardoen, T., 2014. Heterogenous void growth revealed by in situ 3-D X-ray microtomography using automatic cavity tracking. Acta Materialia 63, 130 { 139. Lee, E. H., Liu, D. T., 1967. Finite-strain elastic - plastic theory with appli- cation to plane-wave analysis. Journal of Applied Physics 38, 19{27. Lele, S. P., Anand, L., 2008. A small-deformation strain-gradient theory for isotropic viscoplastic materials. Philosophical Magazine 88, 3655{3689. Li, H., Sun, X., Yang, H., 2016. A three-dimensional cellular automata- crystal plasticity nite element model for predicting the multiscale interac- tion among heterogeneous deformation, drx microstructural evolution and mechanical responses in titanium alloys. International Journal of Plasticity 87, 154 { 180. Lim, H., Abdeljawad, F., Owen, S. J., Hanks, B. W., Foulk, J. W., Battaile, C. C., 2016. Incorporating physically-based microstructures in materials modeling: Bridging phase eld and crystal plasticity frameworks. Mod- elling and Simulation in Materials Science and Engineering 24, 045016. Liu, C., Shanthraj, P., Diehl, M., Roters, F., Dong, S., Dong, J., Ding, W., Raabe, D., 2018. An integrated crystal plasticityphase eld model for spatially resolved twin nucleation, propagation, and growth in hexagonal materials. International Journal of Plasticity 106, 203 { 227. Liu, Y., Li, N., Arul Kumar, M., Pathak, S., Wang, J., McCabe, R. J., Mara, N. A., Tom e, C. N., 2017. Experimentally quantifying critical stresses asso- ciated with basal slip and twinning in magnesium using micropillars. Acta Materialia 135, 411{ 421. Liu, Y., Ponte-Castaneda, ~ P., 2004. Second-order theory for the e ective behavior and eld uctuations of polycrystals. Journal of the Mechanics and Physics of Solids 52, 467{495. LLorca, J., Gonz alez, C., Molina-Aldaregu a, J. M., Segurado, J., Seltzer, R., Sket, F., Rodr guez, M., S adaba, S., Munoz, ~ R., Canal, L. P., 2011. Multiscale modeling of composite materials: a roadmap towards virtual testing. Advanced Materials 23, 5130{5147. 118 LLorca, J., Segurado, J., 2004. Three-dimensional multiparticle cell simula- tions of deformation and damage in sphere-reinforced composites. Materi- als Science and Engineering A 365, 267{274. Lucarini, S., Segurado, J., 2018. On the accuracy of spectral solvers for mi- cromechanics based fatigue modeling. Computational Mechanics. Ma, A., Roters, F., 2004. A constitutive model for fcc single crystals based on dislocation densities and its application to uniaxial compression of alu- minium single crystals. Acta Materialia 52, 3603 { 3612. Ma, A., Roters, F., Raabe, D., 2006. A dislocation density based consti- tutive model for crystal plasticity fem including geometrically necessary dislocations. Acta Materialia 54, 2169{2179. Mahajan, S., Williams, D. F., 1973. Deformation twinning in metals and alloys. International Metallurgical Reviews 18, 43{61. Malyar, N., Micha, J., Dehm, G., Kirchlechner, C., 2017. Size e ect in bi- crystalline micropillars with a penetrable high angle grain boundary. Acta Materialia 129, 312 { 320. Mandal, S., Lao, J., Donegan, S., Rollett, A. D., 2018. Generation of statis- tically representative synthetic three-dimensional microstructures. Scripta Materialia 146, 128 { 132. Manonukul, A., Dunne, F. P. E., 2004. High{ and low{cycle fatigue crack initiation using polycrystal plasticity. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 460, 1881{ Mareau, C., Daymond, M. R., 2016. Micromechanical modelling of twin- ning in polycrystalline materials: Application to magnesium. International Journal of Plasticity 85, 156 { 171. Masson, R., Bornert, M., Suquet, P., Zaoui, A., 2000. Ane formulation for the prediction of the e ective properties of nonlinear composites and polycrystals. Journal of the Mechanics and Physics of Solids 48, 1203{1227. Mayeur, J., Beyerlein, I., Bronkhorst, C., Mourad, H., 2015. Incorporat- ing interface a ected zones into crystal plasticity. International Journal of Plasticity 65, 206 { 225. 119 McColl, I. R., Ding, J., Leen, S. B., 2004. Finite element simulation and experimental validation of fretting wear. Wear 256, 1114 { 1127. McDowell, D. L., Dunne, F. P. E., 2010. Microstructure-sensitive computa- tional modeling of fatigue crack formation. International Journal of Fatigue 32, 1521 { 1542. McDowell, D. L., Gall, K., Horstemeyer, M. F., Fan, J., 2003. Microstructure- based fatigue modeling of cast a356-t6 alloy. Engineering Fracture Mechan- ics 70, 49 { 80. Mecking, H., Kocks, U., 1981. Kinetics of ow and strain-hardening. Acta Metallurgica 29, 1865 { 1875. M eric, L., Poubanne, P., Cailletaud, G., 1991. Single crystal modeling for structural calculations: Part 1 model presentation. ASME Journal of En- gineering Materials and Technology 113, 162{170. Michel, J. C., Moulinec, H., Suquet, P., 1999. E ective properties of com- posite materials with periodic microstructure: a computational approach. Compuational Methods in Applied Mechanics and Engineering 172, 109{ Michel, J. C., Moulinec, H., Suquet, P., 2000. A computational method based on augmented lagrangians and fast Fourier transforms for composites with high contrast. Computer Modeling in Engineering and Sciences 1, 79{88. Michel, J.-C., Suquet, P., Mar 2016. A model-reduction approach to the micromechanical analysis of polycrystalline materials. Computational Me- chanics 57, 483{508. Miehe, C., 1996. Exponential map algorithm for stress updates in anisotropic multiplicative elastoplasticity for single crystals. International Journal for Numerical Methods in Engineering 39, 3367{3390. Miehe, C., Schotte, J., Lambrecht, M., 2002. Homogenization of inelastic solid materials at nite strains based on incremental minimization prin- ciples. application to the texture analysis of polycrystals. Journal of the Mechanics and Physics of Solids 50, 2123 { 2167. 120 Miehe, C., Schroder, J., Schotte, J., 1999. Computational homogenization analysis in nite plasticity simulation of texture development in polycrys- talline materials. Computational Methods in Applied Mechanics and En- gineering 171, 387{418. Mika, D. P., Dawson, P. R., 1999. Polycrystal plasticity modeling of intracrys- talline boundary textures. Acta Materialia 47, 1355 { 1369. Molinari, A., Canova, G. R., Ahzi, S., 1987a. Self-consistent approach of the large deformation polycrystal viscoplasticity. Acta Metallurgica 35, 2983{ Molinari, A., Canova, G. R., Ahzi, S., 1987b. A self consistent approach of the large deformation polycrystal viscoplasticity. Acta Metallurgica 35, 2983{2994. Monchiet, V., Bonnet, G., 2013. Numerical homogenization of nonlinear com- posites with a polarization-based t iterative scheme. Computational Ma- terials Science 79, 276 { 283. Mori, T., Tanaka, K., 1973. Average stress in the matrix and average elastic energy of materials with mis tting inclusions. Acta Metallurgica et Mate- rialia 21, 571{574. Motz, C., Schoeberl, T., Pippan, R., 2005. Mechanical properties of micro- sized copper bending beams machined by the focused ion beam technique. Acta Materialia 53, 4269 { 4279. Moulinec, H., Silva, F., 2014. Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials. International Journal of Numerical Methods in Engineering 97, 960{985. Moulinec, H., Suquet, P., 1994. A fast numerical method for computing the linear and nonlinear mechanical properties of composites. Comptes Rendus de l'Acad emie des Sciences 318, 1417{1423. Moulinec, H., Suquet, P., 1998. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Compua- tional Methods in Applied Mechanics and Engineering 157, 69{94. 121 Mueller, T., Kusne, A., Ramprasad, R., 2017. Machine learning in materials science: Recent progress and emerging applications. Reviews in Computa- tional Chemistry 29, 1866 { 273. Mura, T., 1987. Micromechanics of defects in solids. Martinus-Nijho Pub- lishers. Musinski, W. D., McDowell, D. L., 2012. Microstructure-sensitive probabilis- tic modeling of HCF crack initiation and early crack growth in Ni-base su- peralloy IN100 notched components. International Journal of Fatigue 37, 41 { 53. Musinski, W. D., McDowell, D. L., 2016. Simulating the e ect of grain bound- aries on microstructurally small fatigue crack growth from a focused ion beam notch through a three-dimensional array of grains. Acta Materialia 112, 20 { 39. Nabarro, F. R. N., 1947. Dislocations in a simple cubic lattice. Proceedings of the Physical Society 59, 256. Nayyeri, G., Poole, W. J., Sinclair, C. W., Zae erer, S., 2017. The role of indenter radius on spherical indentation of high purity magnesium loaded nearly parallel to the c-axis. Scripta Materialia 137, 119 { 122. Nemat-Nasser, S., Hori, M., 1993. Micromechanics: Overall Properties of Heterogeneous Materials. Elsevier Science Publishers. Nemcko, M. J., Wilkinson, D. S., 2016. On the damage and fracture of com- mercially pure magnesium using x-ray microtomography. Materials Science and Engineering: A 676, 146 { 155. Neper, 2018. http://neper.sourceforge.net. Nguyen, V.-D., Bchet, E., Geuzaine, C., Noels, L., 2012. Imposing peri- odic boundary condition on arbitrary meshes by polynomial interpolation. Computational Materials Science 55, 390 { 406. Niezgoda, S. R., Turner, D. M., Fullwood, D. T., Kalidindi, S. R., 2010. Optimized structure based representative volume element sets re ecting the ensemble-averaged 2-point statistics. Acta Materialia 58, 4432 { 4445. 122 Niordson, C. F., Kysar, J. W., 2014. Computational strain gradient crystal plasticity. Journal of the Mechanics and Physics of Solids 62, 31 { 47. Nix, W. D., Gao, H., 1998. Indentation size e ects in crystalline materials: A law for strain gradient plasticity. Journal of the Mechanics and Physics of Solids 46, 411 { 425. Nye, J. F., 1953. Some geometrical relations in dislocated crystals. Acta metallurgica 1, 153{162. Ohno, N., Wang, J.-D., 1993. Kinematic hardening rules with critical state of dynamic recovery, part II: Application to experiments of ratcheting be- havior. International Journal of Plasticity 9, 391 { 403. Olson, G. B., 2013. Genomic materials design: The ferrous frontier. Acta Materialia 61, 771{781. Orowan, E., 1934. Zur kristallplastizit at. Zeitschrift fur Physik 89, 605 { 613. Ostien, J., Garikipati, K., 2008. Discontinuous galerkin method for an incompatibility-based strain gradient plasticity theory. In: Reddy, B. D. (Ed.), IUTAM Symposium on Theoretical, Computational and Modelling Aspects of Inelastic Media. Springer, pp. 217{226. Ostoja-Starzewski, M., 2006. Material spatial randomness: From statistical to representative volume element. Probabilistic Engineering Mechanics 21, 112 { 132. Otsuka, T., Brenner, R., Bacroix, B., 2018. Fft-based modelling of transfor- mation plasticity in polycrystalline materials during di usive phase trans- formation. International Journal of Engineering Science 127, 92 { 113. Ozturk, D., Shahba, A., Ghosh, S., 2016. Crystal plasticity fe study of the e ect of thermo-mechanical loading on fatigue crack nucleation in titanium alloys. Fatigue and Fracture of Engineering Materials and Structures 39, 752{769. Panteghihi, A., Bardella, L., 2016. On the nite element implementation of higher-order gradient plasticity, with focus on theories based on plastic distortion incompatibility. Computer Methods in Applied Mechanics and Engineering 310, 810{865. 123 Parton, V. Z., Buryachenko, V. A., 1990. Stress uctuations in elastic com- posites. Soviet Physics Doklady 35, 191{193. Patel, D. K., Kalidindi, S. R., 2017. Estimating the slip resistance from spherical nanoindentation and orientation measurements in polycrystalline samples of cubic metals. International Journal of Plasticity 92, 19{30. Pei, L., Yun-Chang, X., Qing, L., 2011. Plastic anisotropy and fracture be- haviour of AZ31 magnesium alloy. Transactions of Nonferrous Metals So- ciety of China 21, 880{884. Peierls, R., 1940. The size of a dislocation. Proceedings of the Physical Society 52, 34. Peirce, D., 1983. Shear band bifurcations in ductile single crystals. Journal of the Mechanics and Physics of Solids 31, 133 { 153. Peirce, D., Asaro, R. J., Needleman, A., 1982a. An analysis of nonuniform and localized deformation in ductile single crystals. Acta Metallurgica 30, 1087{1119. Peirce, D., Asaro, R. J., Needleman, A., 1982b. An analysis of nonuniform and localized deformation in ductile single crystals. Acta Metallurgica 30, 1087 { 1119. Peirce, D., Asaro, R. J., Needleman, A., 1983. Material rate dependence and localized deformation in crystalline solids. Acta Metallurgica 31, 1951 { Peric, D., de Souza Neto, E. A., Feijo, R. A., Partovi, M., Molina, A. J. C., 2011. On micro-to-macro transitions for multi-scale analysis of non-linear heterogeneous materials: uni ed variational basis and nite element im- plementation. International Journal for Numerical Methods in Engineering 87, 149{170. Petch, N. J., 1953. The cleavage strength of polycrystals. Journal of the Iron and Steel Institute 174, 25{28. Pineau, A., McDowell, D. L., Busso, E. P., Antolovich, S. D., 2016. Failure of metals II: Fatigue. Acta Materialia 107, 484 { 507. 124 Ponte-Castaneda, ~ P., 2002. Second-order homogenization estimates for non- linear composites incorporating eld uctuations: I- Theory. Journal of the Mechanics and Physics of Solids 50, 737 { 757. Ponte-Castaneda, ~ P., 1991. The e ective mechanical properties of non-linear composites. Journal of the Mechanics and Physics of Solids 37, 45{71. Ponte-Castaneda, ~ P., 1996. Exact second-order estimates for the e ective mechanical properties of nonlinear composite materials. Journal of the Me- chanics and Physics of Solids 44, 1757{1788. Proudhon, H., Li, J., Reischig, P., Gueninchault, N., Forest, S., Ludwig, W., 2015. Coupling di raction contrast tomography with the nite element method. Advanced Engineering Materials 18, 903 { 912. Proust, G., Tom e, C., Kaschner, G. C., 2007. Modeling texture, twinning and hardening evolution during deformation of hexagonal materials. Acta Materialia 55, 2137{2148. Przybyla, C., McDowell, D., 2012. Microstructure-sensitive extreme-value probabilities of high-cycle fatigue for surface vs. subsurface crack formation in duplex Ti-6Al-4V. Acta Materialia 60, 293 { 305. Przybyla, C. P., McDowell, D. L., 2010. Microstructure-sensitive extreme value probabilities for high cycle fatigue of Ni-base superalloy IN100. In- ternational Journal of Plasticity 26, 372 { 394. Przybyla, C. P., Musinski, W. D., Castelluccio, G. M., McDowell, D. L., 2013. Microstructure-sensitive HCF and VHCF simulations. International Journal of Fatigue 57, 9 { 27. Pushkareva, M., Adrien, J., Maire, E., Segurado, J., Llorca, J., Weck, A., 2016. Three-dimensional investigation of grain orientation e ects on void growth in commercially pure titanium. Materials Science and Engineering: A 671, 221 { 232. Raabe, D., Sachtleber, M., Zhao, Z., Roters, F., Zae erer, S., 2001. Microme- chanical and macromechanical e ects in grain scale polycrystal plasticity experimentation and simulation. Acta Materialia 49, 3433 { 3441. 125 Radavich, J. F., 1989. The physical metallurgy of cast and wrought alloy 718. In: Loria, E. A. (Ed.), Superalloy 718-Metallurgy and Applications. TMS, pp. 229{240. Redenbach, C., Shklyar, I., Andr a, H., 2012. Laguerre tessellations for elas- tic sti ness simulations of closed foams with strongly varying cell sizes. International Journal of Engineering Science 50, 70 { 78. Reina, C., Conti, S., 2014. Kinematic description of crystal plasticity in the nite kinematic framework: A micromechanical understanding of F=FeFp. Journal of the Mechanics and Physics of Solids 67, 40 { 61. Reina, C., Schlmerkemper, A., Conti, S., 2016. Derivation of F=FeFp as the continuum limit of crystalline slip. Journal of the Mechanics and Physics of Solids 89, 231 { 254. Rice, J., 1971. Inelastic constitutive relations for solids: An internal-variable theory and its application to metal plasticity. Journal of the Mechanics and Physics of Solids 19, 433{455. Richards, A. W., Lebensohn, R. A., Bhattacharya, K., 2013. Interplay of martensitic phase transformation and plastic slip in polycrystals. Acta Ma- terialia 61, 4384 { 4397. Rodr guez-Gal an, D., Sabirov, I., Segurado, J., 2015. Temperature and stain rate e ect on the deformation of nanostructured pure titanium. Interna- tional Journal of Plasticity 70, 191 { 205. Rodr guez-Gal an, D., Segurado, J., Romero, I., 2017. working document. Roters, F., Eisenlohr, P., Bieler, T. R., Raabe, D., 2010a. Crystal Plasticity Finite Element Methods: in Materials Science and Engineering. Wiley- VCH. Roters, F., Eisenlohr, P., Hantcherli, L., Tjahjanto, D. D., Bieler, T. R., Raabe, D., 2010b. Overview of constitutive laws, kinematics, homogeniza- tion and multiscale methods in crystal plasticity nite-element modeling: theory, experiments, applications. Acta Materialia 58, 1152{1211. Sachs, G., 1928. On the derivation of a condition of owing. Verein Deutscher Ingenieure 72, 734 { 736. 126 S anchez-Mart n, R., P erez-Prado, M. T., Segurado, J., Bohlen, J., Guti errez- Urrutia, I., LLorca, J., Molina-Aldareguia, J. M., 2014a. Measuring the critical resolved shear stresses in Mg alloys by instrumented nanoindenta- tion. Acta Materialia 71, 283 { 292. S anchez-Mart n, R., P erez-Prado, M. T., Segurado, J., Bohlen, J., Guti errez- Urrutia, I., LLorca, J., Molina-Aldareguia, J. M., 2014b. Measuring the critical resolved shear stresses in Mg alloys by instrumented nanoindenta- tion. Acta Materialia 71, 283 { 292. Sangid, M. D., 2013. The physics of fatigue crack initiation. International Journal of Fatigue 57, 58 { 72. Sangid, M. D., Maier, H. J., Sehitoglu, H., 2011. An energy-based microstruc- ture model to account for fatigue scatter in polycrystals. Journal of the Mechanics and Physics of Solids 59, 595 { 609. Schmidt, W. L. S., Lauridsen, E. M., Poulsen, H. F., 2008. X-ray di rac- tion contrast tomography: A novel technique for three-dimensional grain mapping of polycrystals. 1. Direct beam case. Journal of Applied Crystal- lography 41, 302 { 309. Segurado, J., Gonz alez, C., LLorca, J., 2003. A numerical investigation of the e ect of particle clustering on the mechanical properties of composites. Acta Materialia 51, 2355{2369. Segurado, J., Lebensohn, R., Llorca, J., Tom e, C., 2012. Multiscale model- ing of plasticity based on embedding the viscoplastic self-consistent for- mulation in implicit nite elements. International Journal of Plasticity 28, 124{140. Segurado, J., LLorca, J., 2002. A numerical approximation to the elastic properties of sphere-reinforced composites. Journal of the Mechanics and Physics of Solids 50, 2107{2121. Segurado, J., LLorca, J., 2005. Computational micromechanics of composites: The e ect of particle spatial distribution. Mechanics of Materials 38, 873{ 127 Segurado, J., LLorca, J., 2013. Simulation of the deformation of polycrys- talline nanostructured Ti by computational homogenization. Computa- tional Materials Science 76, 3{11. Sencer, B. H., Maloy, S. A., Gray, G. T., 2005. The in uence of shock-pulse shape on the structure/property behavior of copper and 316l austenitic stainless steel. Acta Materialia 53, 3293 { 3303. Shahba, A., Ghosh, S., 2016. Crystal plasticity FE modeling of Ti alloys for a range of strain-rates. Part I: A uni ed constitutive model and ow rule. International Journal of Plasticity 87, 48 { 68. Shenoy, M., 2006. Constitutive Modeling and Life Prediction in Ni-Base Su- peralloys. Ph D thesis, Georgia Institute of Technology. Shenoy, M., Zhang, J., McDowell, D., 2007. Estimating fatigue sensitivity to polycrystalline Ni-base superalloy microstructures using a computational approach. Fatigue and Fracture of Engineering Materials and Structures 30, 889{904. Shu, J. Y., Fleck, N. A., 1999. Strain gradient crystal plasticity: size- dependent deformation of bicrystals. Journal of the Mechanics and Physics of Solids 47, 297 { 324. Song, Y. S., Lee, M. R., Kim, J. T., 2005. E ect of grain size for the tensile strength and the low cycle fatigue at elevated temperature of alloy 718 cogged by open die forging press. In: Loria, E. A. (Ed.), Superalloys 718, 625, 706 and Derivatives. TMS, pp. 539{549. Spowart, J. E., 2006. Automated serial sectioning for 3-d analysis of mi- crostructures. Scripta Materialia 55, 5{10. Spowart, J. E., Mullens, H. M., Puchala, B. T., 2003. Collecting and ana- lyzing microstructures in three dimensions: a fully automated approach. JOM 55, 35{37. Staroselsky, A., Anand, L., 2003. A constitutive model for hcp materials deforming by slip and twinning: application to magnesium alloy az31b. International Journal of Plasticity 19, 1843 { 1864. 128 Stelmashenko, N. A., Walls, M. G., Brown, L. M., Milman, Y., 1993. Mi- croindentations on W and Mo oriented single crystals: An STM study. Acta Metallurgica et Materialia 41, 2855 { 2865. Stinville, J., Lenthe, W., Miao, J., Pollock, T., 2016. A combined grain scale elasticplastic criterion for identi cation of fatigue crack initiation sites in a twin containing polycrystalline nickel-base superalloy. Acta Materialia 103, 461 { 473. Stinville, J. C., Vanderesse, N., Bridier, F., Bocher, P., Pollock, T. M., 2015. High resolution mapping of strain localization near twin boundaries in a nickel-based superalloy. Acta Materialia 98, 29 { 42. Stukowski, A., Cereceda, D., Swinburne, T. D., Marian, J., 2015. Thermally- activated non-schmid glide of screw dislocations in w using atomistically- informed kinetic monte carlo simulations. International Journal of Plastic- ity 65, 108 { 130. Su, Y., Zambaldi, C., Mercier, D., Eisenlohr, P., Bieler, T., Crimp, M., 2016. Quantifying deformation processes near grain boundaries in titanium us- ing nanoindentation and crystal plasticity modeling. International Journal of Plasticity 86, 170 { 186. Sulsky, D., Zhou, S. J., Schreyer, H. L., 1995. Application of a particle-in-cell method to solid mechanics. Computational Physics Communications 87, 236{252. Suresh, S., 2012. Fatigue of Materials. Cambridge University Press. Sweeney, C. A., McHugh, P. E., McGarry, J. P., Leen, S. B., 2012. Microme- chanical methodology for fatigue in cardiovascular stents. International Journal of Fatigue 44, 202 { 216. Sweeney, C. A., O'Brien, B., McHugh, P. E., Leen, S. B., 2014a. Experimen- tal characterisation for micromechanical modelling of CoCr stent fatigue. Biomaterials 35, 36 { 48. Sweeney, C. A., OBrien, B., Dunne, F. P. E., McHugh, P. E., Leen, S. B., 2015. Micro-scale testing and micromechanical modelling for high cycle fatigue of CoCr stent material. Journal of the Mechanical Behavior of Biomedical Materials 46, 244 { 260. 129 Sweeney, C. A., OBrien, B., Dunne, F. P. E., McHugh, P. E., Leen, S. B., 2014b. Strain-gradient modelling of grain size e ects on fatigue of cocr alloy. Acta Materialia 78, 341 { 353. Sweeney, C. A., Vorster, W., Leen, S. B., Sakurada, E., McHugh, P. E., Dunne, F. P. E., 2013. The role of elastic anisotropy, length scale and crystallographic slip in fatigue crack nucleation. Journal of the Mechanics and Physics of Solids 61, 1224 { 1240. Tanaka, K., Mura, T., 1981. A dislocation model for fatigue crack initiation. ASME Journal of Applied Mechanics 48, 97{103. Tandon, G., Weng, G., 1988. A theory of particle-reinforced plasticity. Jour- nal of Applied Mechanics 55, 126{135. Taylor, G. I., 1938. Plastic strain in metals. Journal of the Institute of Metals 62, 307 { 324. Taylor, G. I., Elam, C. F., 1923. The distortion of an aluminium crystal dur- ing a tensile test. Proceedings of the Royal Society of London A: Mathe- matical, Physical and Engineering Sciences 102, 643{667. Taylor, G. I., Elam, C. F., 1925. The plastic extension and fracture of alu- minium crystals. Proceedings of the Royal Society of London A: Mathe- matical, Physical and Engineering Sciences 108, 28{51. Thomason, P. F., 1990. Ductile Fracture of Metals. Pergamon Press. Tom e, C., 1999. Self-consistent polycrystal models: a directional compliance criterion to describe grain interactions. Modelling and Simulation in Ma- terial Science and Engineering 7, 723{728. Tom e, C., Canova, G. R., Kocks, U. F., Christodoulou, N., Jonas, J. J., 1984. The relation between macroscopic and microscopic strain hardening in f.c.c. polycrystals. Acta Metallurgica 32 (10), 1637 { 1653. Tom e, C., Lebensohn, R., Kocks, U., 1991. A model for texture development dominated by deformation twinning: Application to zirconium alloys. Acta Metallurgica et Materialia 39 (11), 2667 { 2680. 130 Tom e, C. N., Agnew, S. R., Yoo, M. H., 2001a. Application of texture simu- lation to understanding mechanical behavior of Mg and solid solution alloy containing Li or Y. Acta Materialia 49, 4277{4289. Tom e, C. N., Lebensohn, R. A., , Necker, C. T., 2002. Orientation correla- tions and anisotropy of recrystallized aluminum. Metallurgical and Mate- rials Transactions A 33, 2635 { 2648. Tom e, C. N., Maudlin, P. J., Lebensohn, R. A., Kaschner, G. C., 2001b. Me- chanical response of zirconium: I. Derivation of a polycrystal constitutive law and nite element analysis. Acta Materialia 49, 3085{3096. Torquato, S., 2001. Random heterogeneous materials. Springer. Totry, E., Gonz alez, C., LLorca, J., 2008. Failure locus of ber-reinforced composites under transverse compression and out-of-plane shear. Compos- ites Science and Technology 68, 829{839. Tucker, J. C., Cerrone III, A. R., Ingra ea, A. R., Rollett, A. D., 2015. Crys- tal plasticity nite element analysis for ren88dt statistical volume element generation. Modelling and Simulation in Materials Science and Engineer- ing 23, 035003. Uchic, M. D., Groeber, M., Wheeler, R., Scheltens, F., Dimiduk, D., 2004. Augmenting the 3d characterization capability of the dual beam FIB-SEM. Microscopy and Microanalysis 10, 1136 { 1137. Uchic, M. D., Groeber, M. A., Dimiduk, D. M., Simmons, J. P., 2006. 3d microstructural characterization on nickel superalloys via serial-sectioning using a dual-beam FIB-SEM. Scripta Materialia 55, 23{28. Upadhyay, M. V., Capolungo, L., Taupin, V., Fressengeas, C., Leben- sohn, R. A., 2016. A higher order elasto-viscoplastic model using fast fourier transforms: e ects of lattice curvatures on mechanical response of nanocrystalline metals. International Journal of Plasticity 83, 126{152. Vachhani, S. J., Kalidindi, S. R., 2015. Grain-scale measurement of slip re- sistances in aluminum polycrystals using spherical nanoindentation. Acta Materialia 60, 27{36. 131 Van Houtte, P., Kanjarla, A. K., Van Bael, A., Seefeldt, M., Delannay, L., 2006. Multiscale modelling of the plastic anisotropy and deformation tex- ture of polycrystalline materials. European Journal of Mechanics A/Solids 25, 634{648. Van Houtte, P., Yerra, S. K., Van Bael, A., 2009. The Facet method: A hierarchical multilevel modelling scheme for anisotropic convex plastic po- tentials. International Journal of Plasticity 25, 332{360. Venkatramani, G., Ghosh, S., Mills, M., 2007. A size-dependent crystal plas- ticity nite-element model for creep and load shedding in polycrystalline titanium alloys. Acta Materialia 55, 3971 { 3986. Vondrejc, J., Zeman, J., Marek, I., 2014. An t-based galerkin method for homogenization of periodic media. Computers and Mathematics with Ap- plications 68, 156 { 173. Vondrous, A., Bienger, P., Schreij ag, S., Selzer, M., Schneider, D., Nestler, B., Helm, D., M onig, R., 2015. Combined crystal plasticity and phase- eld method for recrystallization in a process chain of sheet metal production. Computational Mechanics 55, 439{452. Wan, V., MacLachlan, D., Dunne, F., 2016a. Integrated experiment and mod- elling of microstructurally-sensitive crack growth. International Journal of Fatigue 91, 110 { 123. Wan, V. V. C., Jiang, J., MacLachlan, D. W., Dunne, F. P. E., 2016b. Microstructure-sensitive fatigue crack nucleation in a polycrystalline ni superalloy. International Journal of Fatigue 90, 181 { 190. Wang, F., Sandlobes, S., Diehl, M., Sharma, L., Roters, F., Raabe, D., 2014. In situ observation of collective grain-scale mechanics in Mg and Mg-rare earth alloys. Acta Materialia 80, 77 { 93. Wen, W., Borodachenkova, M., Tom e, C. N., Vincze, G., Rauch, E. F., Barlat, F., Gracio, J. J., 2016. Mechanical behavior of Mg subjected to strain path changes: Experiments and modeling. International Journal of Plasticity 73, 171{183. 132 Wenk, H. R., Canova, G. R., Brechet, Y., Flandin, L., 1997. A deformation- based model for recrystallization of anisotropic materials. Acta Materialia 45, 3283 { 3296. Werwer, M., Cornec, A., 2000. Numerical simulation of plastic deformation and fracture in polysynthetically twinned (PST) crystals of TiAl. Compu- tational Materials Science 19, 97{107. Willis, J. R., 1981. Variational and related methods for the overall properties of composites. Advanced Applied Mechanics 21, 1{78. Willot, F., 2015. Fourier-based schemes for computing the mechanical response of composites with accurate local elds. Comptes Rendus M ecanique 34, 232{245. Wu, P. D., Tom e, C. N., Agnew, S. R., Raeisinia, B., Wang, H., 2010. Eval- uation of self-consistent polycrystal plasticity models for magnesium alloy AZ31B sheet. International Journal of Solids and Structures 47, 2905{2917. Xue, X., Righetti, F., Telley, H., Liebling, T. M., Mocellin, A., 1997. The La- guerre model for grain growth in three dimensions. Philosophical Magazine Part B 75, 567{585. Yasi, J. A., Nogaret, T., Trinkle, D. R., Qi, Y., Hector Jr, L. G., Curtin, W. A., 2009. Basal and prism dislocation cores in magnesium: compari- son of rst-principles and embedded-atom-potential methods predictions. Modelling and Simulation in Materials Science and Engineering 17, 055012. Ye mov, S., der Giessen, E. V., Groma, I., 2004. Bending of a single crystal: discrete dislocation and nonlocal crystal plasticity simulations. Modelling and Simulation in Materials Science and Engineering 12, 1069{1086. Yeratapally, S. R., Glavicic, M. G., Hardy, M., Sangid, M. D., 2016. Mi- crostructure based fatigue life prediction framework for polycrystalline nickel-base superalloys with emphasis on the role played by twin bound- aries in crack initiation. Acta Materialia 107, 152 { 167. Zambaldi, C., Yang, Y., Bieler, T. R., Raabe, D., 2012. Orientation informed nanoindentation of ?-titanium: indentation pileup in hexagonal metals deforming by prismatic slip. Journal of Materials Research 27, 356 { 367. 133 Zambaldi, C., Zehnder, C., Raabe, D., 2015. Orientation dependent deforma- tion by slip and twinning in magnesium during single crystal indentation. Acta Materialia 91, 267 { 288. Zecevic, M., Beyerlein, I. J., Knezevic, M., 2017a. Coupling elasto-plastic self-consistent crystal plasticity and implicit nite elements: Applications to compression, cyclic tension-compression, and bending to large strains. International Journal of Plasticity 93, 187 { 211. Zecevic, M., Pantleon, W., Lebensohn, R. A., McCabe, R. J., Knezevic, M., 2017b. Predicting intragranular misorientation distributions in poly- crystalline metals using the viscoplastic self-consistent formulation. Acta Materialia 140, 398{410. Zeman, J., de Geus, T. W. J., Vondrejc, J., Peerlings, R. H. J., Geers, M. G. D., 2017. A nite element perspective on nonlinear FFT-based mi- cromechanical simulations. International Journal for Numerical Methods in Engineering 110, 903{926. Zeman, J., Vond rejc, J., Nov ak, J., Marek, I., 2010. Accelerating a t-based solver for numerical homogenization of periodic media by conjugate gradi- ents. Journal of Computational Physics 229 (21), 8065{8071. Zhang, C., Li, H., Eisenlohr, P., Liu, W., Boehlert, C. J., Crimp, M. A., Bieler, T. R., 2015. E ect of realistic 3d microstructure in crystal plas- ticity nite element analysis of polycrystalline Ti-5Al-2.5Sn. International Journal of Plasticity 69, 21 { 35. Zhang, J., Joshi, S. P., 2012. Phenomenological crystal plasticity modeling an detailed micromechanical investigations of pure magnesium. Journal of the Mechanics and Physics of Solids 60, 945{972.
Journal
Condensed Matter – arXiv (Cornell University)
Published: Apr 7, 2018