Akash Gupta

Mohamed Ben Bettaieb

Farid Abed-Meraim

Surya R Kalidindi

Computationally efficient predictions of crystal plasticity based forming limit diagrams using a spectral database

The present investigation focuses on the development of a fast and robust numerical tool for the prediction of the forming limit diagrams (FLDs) for thin polycrystalline metal sheets using a Taylor-type (full constraints) crystal plasticity model. The incipience of localized necking is numerically determined by the well-known Marciniak-Kuczynski model. The crystal plasticity constitutive equations, on which these computations are based, are known to be highly nonlinear, thus involving computationally very expensive solutions. This presents a major impediment to the wider adoption of crystal plasticity theories in the computation of FLDs. In this work, this limitation is addressed by using a recently developed spectral database approach based on discrete Fourier transforms (DFTs). Significant improvements were made to the prior approach and a new database was created to address this challenge successfully. These extensions are detailed in the present paper. It is shown that the use of the database allows a significant reduction in the computational cost involved in crystal plasticity based FLD predictions (a reduction of about 96% in terms of CPU time).

Introduction

Reliable prediction of the onset of localized necking in thin metal sheets represents a challenging task in the design and manufacturing of sheet metal structural components. During biaxial stretching of sheet metals, the deformation evolves homogeneously at the beginning of the loading. This homogeneous deformation stage is generally followed by diffuse necking, which is characterized by a progressive strain concentration under a smoothly decreasing load. Then, the load abruptly drops with the development of localized necking, which corresponds to the ultimate deformation state that a stretched metal sheet can undergo before failure. For the prediction of localized necking, the well-known forming limit diagram (FLD) concept is most commonly used. This concept was initially introduced by [START_REF] Keeler | Plastic instability and fracture in sheets stretched over rigid punches[END_REF]. Given that the experimental characterization of FLD has proven to be both difficult and expensive, considerable effort has been expended in prior literature [START_REF] Considère | Mémoire sur l'emploi du fer et de l'acier dans les constructions[END_REF][START_REF] Swift | Plastic instability under plane stress[END_REF][START_REF] Hill | On discontinuous plastic states, with special reference to localized necking in thin sheets[END_REF][START_REF] Marciniak | Limit strains in processes of stretch-forming sheet metal[END_REF]Hutchinson and Neale, 1978a) to develop accurate theoretical and/or numerical tools capable of predicting them. Such approaches require the selection of a suitable constitutive framework, for the description of the material state evolution, T together with an instability criterion, for the prediction of the onset of plastic flow localization. These are summarized below:

• Localization criteria: Various analytical as well as numerical criteria have been developed in the literature to predict the onset of plastic strain localization in metals and alloys. These include the imperfection method, initially proposed by [START_REF] Marciniak | Limit strains in processes of stretch-forming sheet metal[END_REF] and subsequently improved by Hutchinson and Neale (1978a), the bifurcation theory developed by Rice and coworkers [START_REF] Rudnicki | Conditions for the localization of deformation in pressure-sensitive dilatant materials[END_REF][START_REF] Stören | Localized necking in thin sheets[END_REF][START_REF] Rice | The localization of plastic deformation[END_REF] and the linear perturbation approach [START_REF] Molinari | Analytical characterization of shear localization in thermo-visco-plastic solids[END_REF][START_REF] Dudzinski | Perturbation analysis of thermoviscoplastic instabilities in biaxial loading[END_REF]. The imperfection approach, which has been extensively used in the literature, is usually referred to as the M-K analysis. It assumes the preexistence of an initial imperfection, which can be of material type (plastic properties) or of geometric nature, in the form of a narrow band across the sheet thickness. This imperfection affects the plastic flow and thus influences the development of localized necking. In the course of loading, the plastic deformation concentrates more and more in the band, ultimately leading to localized thinning in this weakened zone of the sheet. Thanks to its pragmatic character and its extreme flexibility, the initial geometric imperfection approach is widely applied in current literature. This approach is used in this paper.

• Constitutive description: The important role of the material constitutive behavior in the FLD prediction has already been established in prior literature. Accordingly, a wide range of phenomenological constitutive models have been coupled with several localization criteria to investigate the effect on the predicted FLDs of various mechanical features, such as plastic anisotropy [START_REF] Neale | Limit strain predictions for strain-rate sensitive anisotropic sheets[END_REF][START_REF] Cao | Prediction of localized thinning in sheet metal using a general anisotropic yield criterion[END_REF][START_REF] Kuroda | Forming limit diagrams for anisotropic metal sheets with different yield criteria[END_REF][START_REF] Wu | Effect of cube texture on sheet metal formability[END_REF][START_REF] Zhang | Modeling the localized necking in anisotropic sheet metals[END_REF], strain-rate sensitivity [START_REF] Ghosh | Tensile instability and necking in materials with strain hardening and strain-rate hardening[END_REF]Hutchinson and Neale, 1978b;[START_REF] Neale | Limit strain predictions for strain-rate sensitive anisotropic sheets[END_REF][START_REF] Khan | Anisotropic responses, constitutive modeling and the effects of strain-rate and temperature on the formability of an aluminum alloy[END_REF][START_REF] Manopulo | An extended Modified Maximum Force Criterion for the prediction of localized necking under nonproportional loading[END_REF], temperature [START_REF] Khan | Anisotropic responses, constitutive modeling and the effects of strain-rate and temperature on the formability of an aluminum alloy[END_REF][START_REF] Manopulo | An extended Modified Maximum Force Criterion for the prediction of localized necking under nonproportional loading[END_REF] and damage-induced softening [START_REF] Haddag | Strain localization analysis using a large deformation anisotropic elastic-plastic model coupled with damage[END_REF][START_REF] Abed-Meraim | Investigation and comparative analysis of plastic instability criteria: application to forming limit diagrams[END_REF][START_REF] Mansouri | Ductility limit prediction using a GTN damage model coupled with localization bifurcation analysis[END_REF]. However, a major shortcoming in the phenomenological constitutive models is their inability to consider accurately many of the salient microstructural aspects of the material behavior, such as crystallographic texture [START_REF] Hosford | Metal Forming: Mechanics and Metallurgy[END_REF][START_REF] Raabe | Crystal plasticity simulation study on the influence of texture on earing in steel[END_REF]. These limitations represent the main motivation behind the use of micromechanical physics-based constitutive models for the FLD predictions. A salient feature of these physics-based models is their ability to account for the relevant microstructural mechanisms at the single crystal scale, such as slip on specific crystallographic systems, the lattice rotation, and the dislocation motion and multiplication in the description of the strain hardening. The macroscale response is then determined from that of the single crystals via appropriate scale-transition schemes. Such micromechanical modeling approaches have demonstrated remarkable efficiency in predicting the concomitant anisotropic mechanical response in polycrystalline materials and the evolution of the crystallographic texture in finite plastic deformation (e.g., [START_REF] Asaro | Texture development and strain hardening in rate dependent polycrystals[END_REF][START_REF] Bronkhorst | Polycrystalline plasticity and the evolution of crystallographic texture in FCC metals[END_REF][START_REF] Kalidindi | Crystallographic texture evolution in bulk deformation processing of FCC metals[END_REF][START_REF] Beaudoin | A hybrid finite element formulation for polycrystal plasticity with consideration of macrostructural and microstructural linking[END_REF]. Recently, several multiscale models using crystal plasticity constitutive descriptions have been coupled with the M-K analysis in order to predict localized necking and the associated FLDs. It should be noted that the constitutive equations at the single crystal scale can be classified into two main families: rate-dependent (cf. [START_REF] Wu | On crystal plasticity FLD analysis[END_REF][START_REF] Wu | Crystal plasticity forming limit diagram analysis of rolled aluminum sheets[END_REF][START_REF] Wu | Effect of cube texture on sheet metal formability[END_REF][START_REF] Boudeau | Computational prediction of the localized necking in sheet forming based on microstructural material aspects[END_REF][START_REF] Mcginty | Application of multiscale crystal plasticity models to forming limit diagrams[END_REF][START_REF] Inal | Forming limit comparisons for FCC and BCC sheets[END_REF][START_REF] Neil | Crystal plasticity-based forming limit prediction for non-cubic metals: application to Mg alloy AZ31[END_REF][START_REF] Signorelli | Predictions of forming limit diagrams using a rate-dependent polycrystal self-consistent plasticity model[END_REF][START_REF] Lévesque | Numerical modeling of formability of extruded magnesium alloy tubes[END_REF][START_REF] Serenelli | Study of limit strains for FCC and BCC sheet metal using polycrystal plasticity[END_REF][START_REF] Wang | On crystal plasticity formability analysis for magnesium alloy sheets[END_REF][START_REF] Kim | Numerical procedures for predicting localization in sheet metals using crystal plasticity[END_REF][START_REF] Schwindt | Experimental and Visco-Plastic Self-consistent evaluation of forming limit diagrams for anisotropic sheet metals: an efficient and robust implementation of the M-K model[END_REF] and rate-independent approaches (cf. [START_REF] Knockaert | Forming limits predictions using rate-independent polycrystalline plasticity[END_REF][START_REF] Yoshida | Comparison of bifurcation and imperfection analyses of localized necking in rate-independent polycrystalline sheets[END_REF][START_REF] Akpama | Localized necking predictions based on rate-independent self-consistent polycrystal plasticity: bifurcation analysis versus imperfection approach[END_REF]. Elasticity may be considered or disregarded in both approaches. In the current paper, a rigid-viscoplastic framework was adopted to model the mechanical behavior at the single crystal scale. Elasticity is ignored in the constitutive modeling because strain localization occurs at relatively large strains. Furthermore, the deformation follows monotonic strain paths (without elastic unloading) prior to localized necking. Also, in this contribution, attention has been restricted to face-centered cubic (FCC) crystals. The overall mechanical behavior of the polycrystalline aggregate was computed using the Taylor-type scale-transition scheme [START_REF] Taylor | Plastic strains in metals[END_REF]. This particular approach has been shown to provide good predictions for medium to high stacking fault energy cubic materials [START_REF] Asaro | Texture development and strain hardening in rate dependent polycrystals[END_REF][START_REF] Mathur | On modelling the development of crystallographic texture in bulk forming processes[END_REF][START_REF] Bronkhorst | Polycrystalline plasticity and the evolution of crystallographic texture in FCC metals[END_REF]. As explained above, a major advantage of this micromechanical modeling approach compared to the phenomenological one lies in its capability of naturally correlating the microstructure of the polycrystalline metal alloys with their formability. Despite this important benefit, the numerical prediction of FLDs for rigid-viscoplastic materials using the multiscale approach has not been often undertaken for polycrystalline aggregates composed of a large number of single crystals (one typically needs a few thousand single crystals to faithfully represent the polycrystalline aggregate of interest). The main factor limiting the application of the crystal plasticity models in predicting FLDs can be attributed to the considerable computational cost involved. Even when the constitutive equations are developed and implemented efficiently, the time integration of the highly nonlinear single crystal model generally demands high computational resources. Consequently, in recent years, much attention has been devoted to the development of novel techniques that would lead to the same predictions as conventional methods, while consuming much less computational resources. In this field, various numerical strategies have been developed in order to speed up the crystal plasticity computations. These have included the Bunge-Esling approach [START_REF] Bunge | Texture development by plastic deformation[END_REF][START_REF] Li | Evolution of crystal orientation distribution coefficients during plastic deformation[END_REF][START_REF] Kalidindi | Spectral methods for capturing crystallographic texture evolution during large plastic strains in metals[END_REF], direct fast Fourier transform (FFT)-based formulation [START_REF] Lebensohn | N-site modeling of a 3D viscoplastic polycrystal using fast Fourier transform[END_REF][START_REF] Lebensohn | Orientation image-based micromechanical modelling of subgrain texture evolution in polycrystalline copper[END_REF][START_REF] Grennerat | Experimental characterization of the intragranular strain field in columnar ice during transient creep[END_REF] and Fourier-based spectral database approach, which uses generalized spherical harmonics (GSH) for developing efficient spectral representations for functions involved in the crystal plasticity calculations [START_REF] Kalidindi | Spectral calibration of crystal plasticity models[END_REF]. In the current paper, a compact spectral database approach, based on the discrete Fourier transforms (DFTs), is successfully used to speed up the conventional crystal plasticity computations. The viability and the computational benefits of using the DFT approach (for facilitating crystal plasticity solutions) have already been highlighted in prior work (cf. [START_REF] Knezevic | Computationally efficient database and spectral interpolation for fully plastic Taylor-type crystal plasticity calculations of face-centered cubic polycrystals[END_REF][START_REF] Knezevic | Crystal plasticity simulations using discrete Fourier transforms[END_REF][START_REF] Kalidindi | Representation of the orientation distribution function and computation of first-order elastic properties closures using discrete Fourier transforms[END_REF][START_REF] Alharbi | Spectral approaches for the fast computation of yield surfaces and first-order plastic property closures for polycrystalline materials with cubic-Triclinic textures[END_REF][START_REF] Knezevic | A high-performance computational framework for fast crystal plasticity simulations[END_REF][START_REF] Mihaila | Three orders of magnitude improved efficiency with high-performance spectral crystal plasticity on GPU platforms[END_REF][START_REF] Savage | Computer implementations of iterative and non-iterative crystal plasticity solvers on high performance graphics hardware[END_REF][START_REF] Alharbi | Crystal plasticity finite element simulations using a database of discrete Fourier transforms[END_REF][START_REF] Zecevic | Spectral database solutions to elasto-viscoplasticity within finite elements: application to a cobalt-based FCC superalloy[END_REF][START_REF] Eghtesad | Spectral database constitutive representation within a spectral micromechanical solver for computationally efficient polycrystal plasticity modelling[END_REF][START_REF] Knezevic | Crystal plasticity modeling of microstructure evolution and mechanical fields during processing of metals using spectral databases[END_REF]. In particular, the accuracy and computational savings of the DFT database approach compared to the conventional Taylor-type approach have been demonstrated by [START_REF] Knezevic | Crystal plasticity simulations using discrete Fourier transforms[END_REF] and [START_REF] Alharbi | Crystal plasticity finite element simulations using a database of discrete Fourier transforms[END_REF] for both monotonic and non-monotonic strain paths. It was also demonstrated by [START_REF] Alharbi | Spectral approaches for the fast computation of yield surfaces and first-order plastic property closures for polycrystalline materials with cubic-Triclinic textures[END_REF] that the database of DFTs can be used for fast and highly accurate computation of yield surfaces by exploiting the spectral representations of the texture and the stress function. The current paper represents the first application of the DFT spectral database strategy to predict localized necking and the associated FLDs of thin metal sheets, represented as spectral crystal plasticity FLD (SCP-FLD). This new tool is presented and validated in this paper by comparing the predicted SCP-FLDs with the predictions from conventional methods, referred to as CP-FLD (e.g., [START_REF] Wu | On crystal plasticity FLD analysis[END_REF][START_REF] Wu | Crystal plasticity forming limit diagram analysis of rolled aluminum sheets[END_REF][START_REF] Wu | Effect of cube texture on sheet metal formability[END_REF]. We demonstrate that the strategy proposed here leads to major savings of computational cost (96% reduction in the computational cost for FLD predictions).

The rest of the paper will be presented in the following major parts: ✓ The theoretical framework, which is the basis of the developed numerical tool, is presented in Section 2. The rigid-viscoplastic crystal plasticity constitutive framework is first detailed. Thereafter, the equations governing the M-K approach are summarized. ✓ The numerical aspects, related to the resolution and the integration of the constitutive equations developed in Section 2, are detailed in Section 3. Particularly, the details of the conventional method and the DFT approach, required to solve the constitutive equations at the single crystal scale, are expounded. Within the conventional method, the set of equations governing the M-K approach are formulated as a nonlinear mathematical system and solved using the Newton-Raphson iterative scheme. Significant improvements designed and implemented to produce new spectral DFT databases needed for the current work are also detailed in this section. ✓ Comparisons of the predicted FLDs from the new tool (based on the DFT approach) and the conventional approach are presented and discussed in Section 4. ✓ The conclusions of this work are presented in the last section of this paper.

Theoretical framework

Crystal plasticity framework

The mechanical behavior of the sheet is assumed to be described by a rigid-viscoplastic model. Furthermore, the plastic deformation is assumed to be only due to the slip on the crystallographic slip systems. Under these assumptions, the velocity gradient, L, at the single crystal scale can be expressed as

∑ = + = ⊗ * = ⎯ → ⎯⎯⎯ ⎯→ ⎯ γ L W L L m n ; ˙, α N α α α p p 1 s (1)
where * W is the lattice spin tensor, L p is the plastic velocity gradient tensor, γ˙α is the slip rate on the crystallographic system α, N s is the total number of slip systems (equal to 12 for the FCC single crystals), and n The evolution of the slip rate γ˙α is defined by a power law, as follows:

∀ = = = α N γ γ τ s τ τ S P 1,..., : ˙˙sgn ( ) ; : , s α α α m α α α 0 1/ (2)
where S is the deviatoric part of the Cauchy stress tensor σ, while τ α , s α , and P α are the resolved shear stress, the slip resistance, and the symmetric part of the Schmid tensor

⊗ ⎯ → ⎯⎯⎯ ⎯→ ⎯ m n α α
, respectively, associated with slip system α. γ˙0 and m represent the reference value of the slip rates and the strain-rate sensitivity parameter, respectively.

In order to use a vector notation, which is more suitable for an efficient numerical implementation, τ α can be expressed as (equivalent to Eq. ( 2) and noting that S and P α are both traceless)

∀ = = ∼ α N τ S P 1,..., :
. ,

s α α v v (3)
where S v denotes the vector storage of tensor

s α α α α α α α α v 1 1 2 2 1 1 2 2 1 2 2 3 1 3 (4)
The combination of Eqs. sgn ( . ) .

α N α α m α α v 0 1 v v 1/ v v v s (6)
From Eq. ( 6), a condensed matrix (5x5 matrix) form, M, for the microscopic affine modulus can be defined as

∑ = ∂ ∂ = ⊗ ∼ ∼ = - γ m s s M D S P P S P ˙. . α N α α α α α m m v v 0 1 v v v v (1 )/ s (7)
Matrix M will be used in Sections 3.1 and 3.3 instead of the traditional fourth-order tensor form of the affine modulus. This choice greatly improves the computational efficiency of the simulations. The lattice spin tensor * W and the related lattice rotation tensor R are defined as

∑ = = - = * = γ W RR W W W A ˙; ˙, α N α α T p p 1 s ( 8
)
where W is the applied spin tensor (the skew-symmetric part of L), W p is the plastic spin tensor, and A α is the skew-symmetric part of

⊗ ⎯ → ⎯⎯⎯ ⎯→ ⎯ m n α α
. The evolution of the slip resistance s α in terms of the slip rates is described by an isotropic saturation-type law

∑ ⎜ ⎟ ∀ = = ⎛ ⎝ -⎞ ⎠ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ = α N s h s s γ 1,..., : ˙1 ˙, s α α s a β N β 0 1 s (9)
where h 0 , s s and a are hardening parameters.

The full-constraints Taylor model is used to determine the overall mechanical behavior of the polycrystalline aggregate from the mechanical behavior of its microscopic constituents (the single crystals). This model is based on the assumption of a uniform velocity gradient over all the crystals that compose the polycrystalline aggregate. Accordingly, the macroscopic velocity gradient L is assumed to be equal to its microscopic counterparts L (relating to individual grains). By assuming that the microscopic Cauchy stress tensor σ is homogeneous over each grain (or single crystal), the deviatoric part of the macroscopic Cauchy stress tensor S can be related to its microscopic counterpart S i , corresponding to grain i, by

∑ = = v S S, i N i i 1 g (10)
where v i is the volume fraction corresponding to grain i, and N g is the total number of grains that form the polycrystalline aggregate. Analogous to the treatment at the microscopic scale, the condensed form of the macroscopic affine modulus M can be defined as follows:

= ∂ ∂ M D S , v v (11)
where D v denotes the vector storage of the macroscopic strain rate tensor D (the symmetric part of L), and S v denotes the vector storage of S.

Making use of Eqs. ( 7), ( 10) and ( 11), M can be related to the microscopic affine modulus M by the following relation:

∑ = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = - - v M M ( ) . i N i i 1 1 1 g (12)
2.2. Imperfection approach 2.2.1. Generic form for the Cauchy stress and the velocity gradient tensors Consistent with prior literature (Hutchinson et al., 1978), and because the sheets are assumed to be thin, the generalized plane stress assumption is adopted in the M-K analysis. This plane stress condition is generally expressed in terms of the out-of-plane components of the macroscopic Cauchy stress tensor σ by

= = = σ σ σ 0. 13 23 33 (13)
Accordingly, σ can be expressed in the following generic form:

= ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ σ σ σ σ σ 0 0 0 0 0 . 11 12 12 22 (14)
In the current work, we restrict our attention to the prediction of localized necking in polycrystals with orthotropic crystallographic texture. The consideration of orthotropic texture in conjunction with the plane stress state implies that the out-of-plane components of the macroscopic velocity gradient L are equal to zero [START_REF] Knockaert | Rate-independent crystalline and polycrystalline plasticity, application to FCC materials[END_REF][START_REF] Signorelli | Predictions of forming limit diagrams using a rate-dependent polycrystal self-consistent plasticity model[END_REF][START_REF] Yoshida | Comparison of bifurcation and imperfection analyses of localized necking in rate-independent polycrystalline sheets[END_REF][START_REF] Tadano | Plastic flow localization analysis of heterogeneous materials using homogenization-based finite element method[END_REF]

) = = = = L L L L 0.
13 31 23 32

(15)

Furthermore, due to the rigid-viscoplastic behavior considered, the incompressibility condition implies that

=- + L L L (
) .

33 11 22

(16)

By combining Eqs. ( 15) and ( 16), the relevant form of L applicable to our study can be expressed as (17)

= ⎛ ⎝ ⎜ ⎜ - + ⎞ ⎠ ⎟ ⎟ L L L L L L L 0 0 0 0 ( ) .

Formulation of the imperfection approach

The M-K analysis is based on the assumption of a preexisting initial geometric imperfection in the form of a narrow band across the width of the metal sheet, as illustrated in Fig. 1. The width of this band is assumed to be very small and its exact value does not mater, since it is not a relevant parameter in this classical necking modeling approach. Indeed, in this conventional initial imperfection approach, only the relative thickness of the groove is of importance, which is updated all along the computations. Accordingly, the groove shown in Fig. 1 is plotted with a large width only to clearly illustrate the concept of the initial imperfection approach. Quantities outside and inside the band are designated by superscripts H and B, respectively.

For the FLD prediction, the sheet is submitted to in-plane biaxial loading along the axes → x 1 (which can be selected to be aligned with the rolling direction of the sheet) and → x 2 (which can be aligned with the transverse direction of the sheet), as shown in Fig. 1. Furthermore, the sheet is assumed to be very thin, which justifies the adoption of the assumption of plane stress. Hutchinson et al. (1978) have proven the validity of such an assumption in the case of thin sheets. As a consequence of this classical assumption, the band remains normal to the plane of the sheet during the deformation, and the outer surfaces (in the thickness direction) remain in the plane ( → x 1 , → x 2 ). The loading is prescribed so that the strain-path ratio remains constant in the homogeneous zone all along the deformation

= = L L ρ constant. 22 H 11 H (18)
The incompressibility of the sheet implies that the macroscopic velocity gradient in the homogeneous zone has the following form:

= ⎛ ⎝ ⎜ ⎜ -+ ⎞ ⎠ ⎟ ⎟ L ρ L ρ L L 0 0 0 0 0 0 (1 )
. To determine the complete FLD, the strain-path ratio ρ is varied between -1/2 (uniaxial tensile state) and 1 (equibiaxial tensile state).

The M-K approach is also based on the assumption of the uniformity of deformation fields both outside and inside the band, which can be viewed as a conventional idealized modeling approach for the prediction of forming limit diagrams. This assumption has been adopted in all of the contributions in which the initial imperfection approach has been used in conjunction with crystal plasticity multiscale schemes (e.g., [START_REF] Inal | Forming limit comparisons for FCC and BCC sheets[END_REF][START_REF] Knockaert | Forming limits predictions using rate-independent polycrystalline plasticity[END_REF][START_REF] Neil | Crystal plasticity-based forming limit prediction for non-cubic metals: application to Mg alloy AZ31[END_REF][START_REF] Schwindt | Experimental and Visco-Plastic Self-consistent evaluation of forming limit diagrams for anisotropic sheet metals: an efficient and robust implementation of the M-K model[END_REF][START_REF] Signorelli | Predictions of forming limit diagrams using a rate-dependent polycrystal self-consistent plasticity model[END_REF]. We recall that the initial imperfection approach is used to predict the incipience of localized necking, which is material type instability and, as such, independent of the geometric characteristics (in contrast to geometric or structural type instabilities). Consequently, all geometric or structural effects are purposely excluded when this approach is used, as well as the associated structural type instability phenomena (e.g., wrinkling, plastic buckling). In this approach, the band and the safe zone are restricted to remain homogeneous and, to this end, a very small representative volume element (RVE) or a material point is taken to model the mechanical behavior of these zones. At the macroscopic scale, this RVE is viewed as a homogeneous medium where the macroscopic deformation field remains uniform during the loading (not the same for both zones). However, the behavior is strongly heterogeneous at the microscopic level. Based on this consideration, the current thickness e H (resp. e B ) is related to the initial thickness e 0 H (resp. e 0 B ) by the following relation:

= = e e ε e e ε exp ( ) ; exp ( ), H 0 H 33 H B 0 B 33 B ( 20
)
where ε 33 H and ε 33 B are the logarithmic strains outside and inside the band, respectively, along the direction of the normal to the sheet.

The above-mentioned quantities can be used to quantify the initial and current geometric imperfections through the computation of the ratios

= = f e e f e e ;
.

0 0 B 0 H B H (21)
Using Eqs. ( 20) and ( 21), the current imperfection ratio f can be expressed in terms of ε 33 H and ε 33 B as follows:

= - f f ε ε exp ( ). 0 33 B 33 H ( 22
)
The initial (resp. current) imperfection is also characterized by the orientation ψ 0 (resp. ψ) of the unit vector normal to the band, ⎯ → ⎯⎯ N 0 (resp. ⎯ → ⎯⎯ N ), with respect to the rolling direction of the sheet. Both ⎯ → ⎯⎯ N 0 and

⎯ → ⎯⎯

N are constrained to lie in the → x 1 -→ x 2 plane and are equal to ψ ψ (cos , sin ) 0 0 and ψ ψ (cos , sin ), respectively. The current orientation ψ of the band can be expressed in terms of its initial value ψ 0 , the logarithmic strain component ε 11 H , and the strain-path ratio ρ equal to ε ε /

= - ψ ψ ρ ε arctan[tan( ) exp(1 )
].

0 11 H

In addition to Eqs. ( 20)-( 23), the M-K analysis is based on three other sets of equations:

• The kinematic compatibility condition between the band and the homogeneous zone (this condition is restricted to the in-plane components of the macroscopic velocity gradient)

= + ⊗ ⎯→ ⎯ ⎯ → ⎯⎯ L L C N ˙, PS B PS H (24)
where the notation • PS denotes the in-plane part, defined as

⎜ ⎟ ⎛ ⎝ ⎞ ⎠ • • • • 11 12 21 22
, of tensor • (where • in Eq. ( 24) corresponds to tensors L B and L H ). In Eq. ( 24),

⎯→ ⎯

C ˙is a 2D vector representing the jump in the velocity gradient across the boundary between the band and the homogeneous zone.

The fully three-dimensional form of the velocity gradient in the band is complemented starting from Eq. ( 24) by using the incompressibility assumption, which leads to the following expression for tensor L B :

= ⎛ ⎝ ⎜ ⎜ + + -+ - - ⎞ ⎠ ⎟ ⎟ L C N C N C N ρL C N ρ L C N C N L ˙˙0 ˙˙0 0 0 ( 1 ) ˙˙. B 11 H 1 1 1 2 2 1 11 H 2 2 11 H 1 1 2 2
(25)

• The equilibrium equation across the interface between the band and the homogenous zone

→ = → e e σ N σ N. . . B PS B H PS H (26)
• The constitutive equations developed in Section 2.1 for the modeling of the mechanical behavior of polycrystalline aggregates representing both the band and the homogeneous zone.

For practical reasons related to the numerical integration, the equilibrium equation ( 26) is written in the following compact matrix form:

+ = f B σ b 0, . (1/ ) r B ( 27
)
where f is the current imperfection ratio (see Eq. ( 21)), while B, σ r B and b are defined as is equal to t Δ ), by using an implicit/explicit algorithm: implicit for the computation of the deviatoric stress components and explicit for the update of the lattice rotation tensor R as well as the slip resistances s α . In this process, the values of the mechanical variables σ, R and s α (for = α N 1,..., s ) at t n , and those of the material parameters (hardening parameters) are assumed to be known. The aim of this algorithm is to compute the values of σ, R and s α at time + t n 1 . The microscopic velocity gradient L is also assumed to be known and constant over I Δ . The proposed algorithm is very similar in spirit to that developed by [START_REF] Kalidindi | Crystallographic texture evolution in bulk deformation processing of FCC metals[END_REF]. It is defined by the following steps: (i) Compute t P ( ) α n by using the following transformation:

= ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ - + - + ⎤ ⎦ ⎥ N N N N σ σ σ N σ N σ N σ N σ B σ b 0 0 ; ; ( ) ( ) .
∀ = = α N t t t P R P R 1,..., : ( ) ( ). . ( ), s α n n α n 0 T (29)
where the expression of

P α 0 is given in Appendix A. Once t P ( ) α n is computed, vectors t P ( ) α n v and ∼ t P ( ) α n v
can be directly deduced.

(ii) Solve the following set of scalar equations for the five independent components of

+ t S (
)

n v 1 ∑ = - = ∼ ∼ + = + + t γ t t s t t t t T D S P S P P 0 ( ) ˙( ). ( ) ( ) sgn ( ( ). ( )) ( ) , n α N n α n α n m n α n α n 1 v 0 1 v 1 v 1/ v 1 v v s ( 30
)
by using the Newton-Raphson iterative algorithm as follows:

= + + + + - + + t t t t S S M T ( ) ( ) ( ). ( ), k n k n k n k n v ( 1) 1 v ( ) 1 1( ) 1 ( ) 1 (31)
where • k ( ) is the k th Newton-Raphson iteration of • . Note that, for moderate to large values of the strain-rate sensitivity parameter m, the Newton-Raphson procedure seems to be the most efficient method for the solution of the crystal plasticity constitutive equations. Indeed, as shown in Eq. ( 31), the Jacobian matrix used in the application of the Newton-Raphson procedure is equal to the inverse of the microscopic affine modulus determined in Eq. ( 7).

However, for small values of the strain-rate sensitivity parameter m, the set of equations T becomes very stiff, and it may turn out necessary to apply a constraint to the Newton-Raphson corrections (31), which involves a modification in the iterative procedure. The aim of this modification is to improve the convergence of the iterative process as explained below.

Let

+ + γ t ˙( ) α k n ( 1)
1 denote the slip rate of a given slip system α, computed at iteration + k 1 as follows:

∀ = = ∼ ∼ + + + + + + α N γ t γ t t
s α k n k n α n α n m k n α n ( 1) 1 0 v ( 1) 1 v 1/ v ( 1) 1 v (32) If < + + - γ t η ˙( ) s α k n ( 1) 1
1 for all the slip systems α (where η is a threshold value taken as 2 in this work), then the Newton-Raphson correction previously applied in Eq. ( 31) is sufficient and no further modification is required.

Else, correction (31) must be modified by introducing a scalar λ 1 such that:

= + + + + - + + t t λ t t S S M T ( ) ( ) ( ). ( ). k n k n k n k n v ( 1) 1 v ( ) 1 1 1( ) 1 ( ) 1 (33)
The introduction of scalar λ 1 aims to implement a line search strategy to ensure and accelerate the convergence of the global Newton-Raphson method used to solve Eq. ( 30). Indeed, scalar λ 1 should be determined in such a way that

+ + γ t ˙( ) α k n ( 1) 1 , computed by using + + t S ( ) k n v ( 1)
1 as determined from Eq. ( 33), does not exceed the threshold value η for all the slip systems. It must be noted that this threshold value does not have any effect on the value of the slip rates obtained at the convergence of the Newton-Raphson scheme. The initial guess for

+ S t ( ) n v 1 used to start the iterative procedure is taken equal to t S ( ) n v
. Iterations are carried out until the L 2 norm of

+ + t T ( ) k n ( 1)
1 becomes less than 10 -8 s -1 .

(iii) Store the microscopic affine modulus M:

∑ = ⊗ ∼ ∼ + = + - t γ m t t s t t t s t M P P S P ( ) ˙( ) ( ) ( ) ( ). ( ) ( ) . n α N α n α n α n n α n α n m m 1 0 1 v v v 1 v (1 )/ s ( 34
)
(iv) Update the lattice rotation tensor R:

= - + t t t R W W R ( ) exp(Δ ( )). ( ), n n 1 p ( 35
)
where W is the skew-symmetric part of L, and W p is the plastic part of W, which is determined as follows:

= ∑ ∀ = = = + γ t t α N t t t W A A R A R ˙( ) ( ) 1,..., : ( ) ( ). . ( ). α N α n α n s α n n n p 1 1 0 α T s (36)
the expression of tensor A α 0 is given in Appendix A. Furthermore, it is noteworthy that the slip rates

+ γ t ˙( )
α n 1 for the different slip systems are computed at the convergence of the iterative procedure of Step (ii).

(v) Compute the new values of the slip resistances

+ s t ( ) α n 1 : ⎜ ⎟ ⎛ ⎝ ∑ ⎞ ⎠ ⎜ ⎟ ∀ = = + ⎛ ⎝ - ⎞ ⎠ + = + α N s t s t h s t s t γ t 1,..., : ( ) ( ) 1 ( ) Δ ˙( ) . s α n α n α n s a β N β n 1 0 1 1 s ( 37
)
From Steps (iv) and (v), it can be seen that the lattice rotation R and the slip resistance s α are obtained by an explicit integration scheme. It was found that the results obtained with explicit update of R and s α were not significantly different from those obtained by a fully implicit integration. Consequently, the gain in terms of accuracy allowed by a fully implicit scheme is minimal, while the latter requires a much larger computation time. Indeed, this fully implicit scheme is generally based on two-level methods, as has been detailed in [START_REF] Kalidindi | Crystallographic texture evolution in bulk deformation processing of FCC metals[END_REF], with an external level to update the lattice rotation and the slip resistance, and an internal level to compute the deviatoric stress.

After applying the previous integration algorithm to all the single crystals that compose the polycrystalline aggregate, the deviatoric part of the macroscopic Cauchy stress tensor S should be calculated from its microscopic counterparts S i by using Eq. ( 10). Thereafter, the macroscopic Cauchy stress tensor σ can be determined, under the plane stress condition, through the following relation:

= -S σ S I . 33 2 (38)
When the plane stress assumption is adopted (which is the case in Section 3.3), the 2D macroscopic Cauchy stress tensor σ PS can be obtained from the components of the 3D Cauchy stress tensor σ by the following relation:

= ⎛ ⎝ ⎞ ⎠ σ σ σ σ σ . PS 11 12 12 22 (39)
Likewise, the macroscopic affine modulus M can be obtained from its microscopic counterparts M i , determined in the third step of the previous algorithm, using Eq. ( 12). This modulus M will be used in Section 3.3.2 to integrate the equations governing the imperfection approach.

Crystal plasticity computations using database of DFTs

Spectral analysis is the process of identifying embedded frequencies in data. For discrete data, the computational basis of spectral analysis is discrete Fourier transforms (DFTs). DFTs transform time-or space-based data into frequency-based data [START_REF] Adams | Microstructure-sensitive Design for Performance Optimization[END_REF]. Recently, a new DFT database based spectral approach has been developed [START_REF] Knezevic | Crystal plasticity simulations using discrete Fourier transforms[END_REF][START_REF] Alharbi | Crystal plasticity finite element simulations using a database of discrete Fourier transforms[END_REF][START_REF] Knezevic | Crystal plasticity modeling of microstructure evolution and mechanical fields during processing of metals using spectral databases[END_REF] to solve the crystal plasticity constitutive equations, which reduces the computational cost of performing crystal plasticity calculations tremendously (up to two orders of magnitude) by allowing compact representation and fast retrieval of crystal plasticity solutions. This approach allows one to reconstruct directly the solutions for the main functions of the conventional crystal plasticity theory, for any given crystal orientation and applied deformation mode, using only a small number of DFT terms (called dominant DFTs). Spectral representations are established for the following three main functions required for crystal plasticity solutions: (i) the deviatoric stress tensor g S L ( , ), (ii) the plastic spin tensor g W L ( , ) p , and (iii) the total shear rate

∑ = γ g L ˙( , ) α N α 1 s
, where g denotes the crystal lattice orientation represented using Bunge-Euler angle and L is the applied velocity gradient tensor. We can determine the lattice spin tensor W* from the plastic spin tensor W p using Eq. ( 8), and the Cauchy stress tensor σ to be used for FLD prediction from the deviatoric stress tensor S under the plane stress condition using ( = -S σ S I 33 2 ). The domain of the three main functions to be obtained is the product space of the orientation space and the deformation mode space. Deformation mode space includes all possible L vectors and can be parameterized as (Van Houtte, 1994): can be effectively converted to g p , where g p is the crystal orientation with respect to the e { } j p reference frame [START_REF] Knezevic | Crystal plasticity simulations using discrete Fourier transforms[END_REF]. As a result of the above parametrization, the functions of interest can then be expressed in terms of ( e { } j p , θ, ε˙) instead of L. For computing DFTs, we need to define periodic domains of all independent variables, orientation space (forg p ) and the deformation mode space (for θ). Periodic domain of orientation space g p is identified using Bunge-Euler angles as

= + = =∑ ⊗ = - = + = - = ( ) ( ) ε ε D D θ D θ D θ L D W
∈ ∈ ∈ φ Φ φ ( [0, 2π), [0, 2π), [0, 2π)) 1 2
. For θ, from Fig. 2 below we can identify periodic domain as ∈ θ [0, 2π/3). These periodic domains are used for computing DFTs on a uniform 4-D grid of φ Φ φ θ , , ,

.

The spectral databases for the functions of interest described above are computed on a uniform grid using a model based on conventional crystal plasticity theory described earlier and stored in the form of Fourier coefficients. Interpolation techniques can be used to recover the function values at any other location of interest (i.e., not on a grid point).

To use the DFT spectral database (SCP) model to predict FLDs in this work, the strategies developed in prior work [START_REF] Knezevic | Crystal plasticity simulations using discrete Fourier transforms[END_REF][START_REF] Alharbi | Crystal plasticity finite element simulations using a database of discrete Fourier transforms[END_REF] were significantly improved. This was needed as the DFT databases developed earlier did not exhibit the accuracy needed for the present study. The following specific improvements were implemented to produce new spectral databases for this work:

(1) In new databases, the discretization of a uniform grid in the product space defined by the Bunge-Euler space (φ Φ φ , ,

2 ) and the deformation mode space (θ) was reduced to one degree instead of three degrees used in the old databases. The new databases produced with the one-degree grid provided significantly more accurate results. Improvements on using one-degree grid compared to three-degree grid come from the fact that a finer sampling grid provides more accurate representation of the function in the DFT space. Although there is an additional computation cost in calculating the database on the finer one-degree grid, it should be remembered that it is only a one-time expense. The cost of recovering the solutions from this new database is roughly the same as the earlier database, which is still about two orders of magnitude faster than the conventional computations. The number of simulations required for generating DFT database was reduced by taking advantage of the symmetry relations for cubic crystals [START_REF] Adams | Microstructure-sensitive Design for Performance Optimization[END_REF] and the mirror symmetry evident within the periodic domain of θ (see Fig. 2).

(2) The spectral databases were generated for all six components of the deviatoric stress, unlike the earlier database that addressed only the five independent components. More specifically, in the old database there was no separate database for S 33 and it was computed using relation

=- + S S S ( ) 33 11
22 to ensure that the trace of deviatoric stress tensor is zero. In this work, it was discovered that this strategy resulted in a larger error for S 33 , as a consequence of the addition of the truncation errors for S 11 and S 22 . This problem was addressed by creating a separate database for S 33 . For imposing the condition that the trace of the deviatoric stress tensor has to be zero, the hydrostatic part was subtracted from each of the normal stress components. After this correction, we obtained much more accurate predictions of the entire deviatoric stress tensor from spectral database as will be evident from results below.

(3) In the recovery of results from the database for new values of inputs (combinations of grain orientations and the imposed velocity gradient tensor), trigonometric interpolation was used instead of the spectral interpolation used in prior work [START_REF] Knezevic | Crystal plasticity simulations using discrete Fourier transforms[END_REF][START_REF] Alharbi | Crystal plasticity finite element simulations using a database of discrete Fourier transforms[END_REF]. Spectral interpolation introduced certain rounding errors, which were eliminated by the use of trigonometric interpolation, which is a continuous interpolation scheme [START_REF] Zecevic | Spectral database solutions to elasto-viscoplasticity within finite elements: application to a cobalt-based FCC superalloy[END_REF]. This change resulted in significant improvement in the accuracy of the crystal plasticity computations.

The main functions of crystal plasticity solutions are recovered as continuous functions of orientation and deformation mode from spectral DFT database with trigonometric interpolation using Inverse Discrete Fourier Transforms (IDFTs) equations in principal frame of D, given by [START_REF] Zecevic | Spectral database solutions to elasto-viscoplasticity within finite elements: application to a cobalt-based FCC superalloy[END_REF]:

, ) ˙, p N N N N α β γ δ αβγδ p m N N N N α β γ δ αβγδ α N α p N N N N α β γ δ αβγδ p 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 φ Φ φ θ Nφ Nφ N Φ N Φ Nφ Nφ N θ N θ iαφ Lφ iβΦ L Φ iγφ Lφ iδθ L θ φ Φ φ θ Nφ Nφ N Φ N Φ Nφ Nφ N θ N θ iαφ Lφ iβΦ L Φ iγφ Lφ iδθ L θ s φ Φ φ θ Nφ Nφ N Φ N Φ Nφ Nφ N θ N θ iαφ Lφ iβΦ L Φ iγφ Lφ iδθ L θ 1 2 1 2 1 2 2 2 2 2 2 2 2 2 2 π 1 1 2 π 2 π 2 2 2 π 1 2 1 2 1 2 2 2 2 2 2 2 2 2 2 π 1 1 2 π 2 π 2 2 2 π 1 2 1 2 1 2 2 2 2 2 2 2 2 2 2 π 1 1 2 π 2 π 2 2 2 π (41) whereN N N N , , , φ Φ φ θ 1 2 (
are the total number of grid points and L L L L , , ,

φ Φ φ θ 1 2
represent the domain of variables in the orientation space (g p or (φ Φ, φ , 1 2 ) and the deformation mode space (θ). For one-degree grid and periodic domain considered here

= = = N N N 360 φ Φ φ 1 2 , = N 120 θ and = = = L L L 360 φ Φ φ 1 2 , = L 120 θ since periodic domain of variables considered are ∈ ∈ ∈ φ Φ φ [0, 2π), [0, 2π), [0, 2π), 1 2
∈ θ [0, 2π/3). When retrieving Fourier coefficients (B , αβγδ G C , αβγδ αβγδ ) referred as DFTs for use in above equations to reconstruct values of above functions, we use only dominant DFTs and also exploit the fact that about half of the DFTs are complex conjugates of other half for all real valued functions. Use of small number of dominant DFTs ensures that even with using one degree discretization of uniform grid in new database, the speed up of crystal plasticity calculations obtained continues to be about two orders of magnitude [START_REF] Knezevic | Crystal plasticity simulations using discrete Fourier transforms[END_REF].

One of the salient aspects of the DFT approach is that the Fourier coefficients (B , αβγδ G C , αβγδ αβγδ ) do not exhibit any path dependence (within the assumption of equal hardening of all slip systems at any selected material point in a crystalline domain) as discussed in prior work [START_REF] Knezevic | Computationally efficient database and spectral interpolation for fully plastic Taylor-type crystal plasticity calculations of face-centered cubic polycrystals[END_REF][START_REF] Knezevic | Crystal plasticity simulations using discrete Fourier transforms[END_REF]. Since the DFT database is applicable to all textures and all deformation histories, the path independence of Fourier coefficients is of paramount importance. Building spectral DFT database is a computationally expensive task, but it is only a one time cost, since the same database can be used for a broad range of crystal plasticity computations on FCC materials deforming by ( 111) [110] slip systems at low homologous temperatures irrespective of the specific details of the hardening laws (as long as one imposes equal hardening on all slip systems). After obtaining the simulation data from conventional crystal plasticity model, it took a computational time of ∼70 h on a 2 GHz clock speed machine with a memory (RAM) of ∼128 GB to generate the complete DFT database for all output variables of interest (a total of ten).

Fig. 3 and Fig. 4 provide example results of deformed texture (pole figures) and deviatoric stress-strain curves comparing the accuracy of rigid-viscoplastic crystal plasticity calculations obtained using new and old DFT databases against conventional CP Taylor-type model for deformation modes corresponding to two extremes of the FLD: uniaxial tension (Fig. 3) and equibiaxial tension (Fig. 4). For this study, the hardening parameters were taken from [START_REF] Kalidindi | Crystallographic texture evolution in bulk deformation processing of FCC metals[END_REF] and = m 0.01. Several studies have been carried out previously to evaluate the minimum number of single crystals that should be used to ensure the representativeness of the polycrystalline aggregate. These studies reveal that the obtained results depend on the boundary conditions applied on the polycrystalline aggregates (e.g., periodic boundary conditions) as well as on the degree of anisotropy of the studied single crystals. In most of the previous literature (for instance, El Houdaigui et al., 2005;[START_REF] Szyndler | Effect of number of grains and boundary conditions on digital material representation deformation under plane strain[END_REF], the number of grains required to ensure the representativeness of the volume element does not exceed 1000. Accordingly, we have chosen a polycrystal with 1000 grains, considering that this number is sufficient to generate a volume element representative of the studied sheets. Indeed, from a variety of numerical experiments, we have observed that beyond 1000 grains, the response of the polycrystalline aggregate, in terms of stress-strain curves, remains almost unchanged. For SCP model, 500 dominant DFTs were chosen for deviatoric stress components (C αβγδ ) and total shear rate (G αβγδ ), and 2000 dominant DFTs for plastic spin components (B αβγδ ) with a strain step of 0.06%. Examining the predictions in both Figs. 3 and4 from the old and new spectral databases, it should be clear that the old database produces significant errors for some of the stress components. In this work, it was observed that this small error in the individual stress components translates to a very large error in the FLD predictions. However, a significant improvement in accuracy can be clearly observed from Figs. 3 and4 when the new DFT database is used instead of the old database. This increased accuracy was found to be essential for meaningful and practically useful predictions of forming limit diagrams (FLDs).

Numerical aspects related to the FLD prediction

Algorithm for the FLD prediction

To predict the FLDs of sheet metals, the following algorithm is applied:

• For each strain-path ratio ρ varying from -1/2 to 1 ( ρ Δ is typically set to 0.1): o For each time step = + t t I [ , ] n n Δ 1 , use the algorithm of Section 3.1 (when the conventional integration scheme is used), or Section 3.2 (when the DFT approach is employed), to integrate the constitutive equations of the polycrystal in the homogeneous zone. For this purpose, the 3D macroscopic velocity gradient L H given by Eq. ( 19) should be used as input for this time integration scheme, with the component L 11 H fixed to one. Consequently, the components L 22 H and L 33 H are equal to ρ and -+ ρ (1 ), re- spectively. After applying the algorithm of Section 3.1, the 2D macroscopic stress σ PS H computed at + t n 1 should be stored. The other out-of-plane components of σ H (σ 13 H , σ 23 H and σ 33 H ) are equal to zero as a consequence of the plane stress conditions (see Eq. ( 14)). Tensor σ PS H will be used in Section 3.3.2 to determine the in-plane velocity gradient in the band L PS B , which fulfills both the compatibility and the equilibrium conditions at the interface between the homogeneous zone and the band. The computation in the homogeneous zone is stopped when (

∫ = = ε L d t t t 11 H 0 11 H
) reaches a value of one. o For each initial band inclination ψ 0 varying from 0°to 90°(with ψ Δ 0 set to 1°or 1.5°). ❖ For each time step

+ t t [ , ]
n n 1 (for convenience, the time steps used to integrate the mechanical behavior in the band are chosen to be the same as those taken for the integration of the mechanical behavior in the homogeneous zone):

• Use the value of σ PS H at + t n 1 (determined in the previous step) to integrate the equations of the M-K analysis (see Section 2.2.2). For the sake of clarity, the steps followed in this integration will be detailed in a separate section (Section 3.3.2).

The incremental integration scheme of Section 3.3.2 is applied until the following condition is met: The strain component ε 11 H , obtained once the criterion ( 42) is satisfied, is considered as the critical strain ε 11 * , which corresponds to the current band inclination ψ and strain-path ratio ρ.

≥ L L / 10. 33 B 33 H (42)
The smallest critical strain ε 11 * predicted by the above algorithm, over all initial angles ψ 0 , and the associated current angle ψ correspond, respectively, to the necking limit strain ε 11 L and the necking band inclination ψ for the current strain-path ratio ρ.

Numerical integration of the equations governing the imperfection approach

The equations corresponding to the M-K approach, detailed in Section 2.2.2, have been integrated in order to determine the velocity gradient in the band L B over the loading increment I Δ , and therefore to check whether the criterion ( 42) is satisfied. To simplify the following numerical developments, the argument + t n 1 is dropped in this section, with the implied understanding that all variables are evaluated at + t n 1 , unless otherwise specified. As explained in the algorithm of Section 3.3.1, tensor σ PS H is known a priori. We briefly recall the equations governing the M-K analysis to be solved in this section:

• The compatibility condition given by Eq. ( 25).

• The equilibrium equation ( 27).

• The constitutive equations of the polycrystalline aggregate in the band.

The above set of equations will be solved by means of the Newton-Raphson iterative method. As L 11 H is set to 1 (see Section 3.3.1), the macroscopic velocity gradient in the band L B , defined by Eq. ( 25), can be re-expressed as follows:

= ⎛ ⎝ ⎜ ⎜ + + -+ - - ⎞ ⎠ ⎟ ⎟ C N C N C N ρ C N ρ C N C N L 1 ˙˙0 ˙˙0 0 0 ( 1 ) ˙˙. B 1 1 1 2 2 1 2 2 1 1 2 2 (43)
On the other hand, as ε 11 H is known and actually equal to + t n 1 (since L 11 H is fixed to 1), the current band orientation ψ (at + t n 1 ) is also known, as defined by Eq. ( 23). Consequently, the components N 1 and N 2 of

⎯ → ⎯⎯

N are known at + t n 1 . The computation of L B is thus reduced to the computation of the components C ˙1 and C ˙2 of the jump vector. In turn, C ˙1 and C ˙2 can be determined by solving the equilibrium equation ( 27). The latter may be viewed as a system of two nonlinear equations, which is function of the unknown vector

→ = C C C ˙( ˙, ˙) 1 2 = = + = ⌢ ⎯→ ⎯ f E E C B σ b 0. ( ˙) . (1/ ) r B (44)
The classical Newton-Raphson method is then used to solve Eq. ( 44). In this process, and for the computation of the analytical expression of the Jacobian matrix ∂ ∂

E C

/ ˙can be obtained by the following chain rule:

∂ ∂ = ∂ ∂ ∂ ∂ ⎯→ ⎯ ⎯ → ⎯ E C E D D C ˙• ˙. r B r B (46) The partial derivative ∂ ∂ ⎯→ ⎯ D C / ṙB
, introduced in Eq. ( 46), can be computed using Eq. ( 43)

∂ ∂ = ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎯→ ⎯ N N N N D C ˙0 0 /2 /2 . r B 1 2 2 1 (47)
Using Eq. ( 44), the matrix ∂ ∂ E D / r B , introduced in Eq. ( 46), can be expressed as 12)), as follows:

∂ ∂ = ∂ ∂ + ⊗ ∂ ∂ f E D B σ D b D • (1/ ) ,
∂ ∂ = ⎡ ⎣ ⎢ ⎢ ⎢ + + + + + + ⎤ ⎦ ⎥ ⎥ ⎥ - - - - - - - - - - - - - - - M M M M M M M M M M M M M M M σ D 2 2 2 2 2 2 .
As to the imperfection ratio f , it can be determined from Eq. ( 22) as follows:

= + - f f ε t tL ε exp( ( ) Δ ), n 0 33 B 33 B 33 H (50)
where ε t ( )

n 33 B
is the 33 component of the logarithmic strain tensor in the band at t n , and ε 33 H is the 33 component of the logarithmic strain tensor in the homogeneous zone at + t n 1 . As D 33 H is constant and equal to -+ρ (1 ), ε 33 H is simply equal to -+

+ ρ t (1 ) n 1 . Finally, ∂ ∂ f D (1/ )/ r
B can be expressed as follows:

∂ ∂ = ⎡ ⎣ ⎤ ⎦ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ f D (1/ ) 0 . f D D D f D D D r B (1/ ) (1/ ) 33 B 33 B 11 B 33 B 33 B 22 B (51)
Using Eq. ( 50),

∂ ∂ f D (1/ )/ 33 B reads ∂ ∂ =- f D t f (1/ ) Δ . 33 B (52)
The incompressibility condition leads to the following relations:

∂ ∂ = - ∂ ∂ = - D D D D 1 ; 1. 33 B 11 B 33 B 22 B (53)
The insertion of Eqs. ( 52) and ( 53) into Eq. ( 51) gives

∂ ∂ = ⎡ ⎣ ⎤ ⎦ f D (1/ ) 0 . t f t f r B Δ Δ (54)
Substituting Eqs. ( 47)-( 49) and ( 54) into Eq. ( 46) gives the analytical Jacobian matrix

∂ ∂ ⎯→ ⎯ E C / ˙.
The Jacobian matrix ∂ ∂ ⎯→ ⎯

E C

/ ˙is analytically determined by using the above developments when the CP-FLD model is used. However, this matrix is evaluated numerically in the case of SCP-FLD model. This latter choice is adopted because the strain-rate sensitivity parameter m is set to small values (typically about 0.01), which leads to a very stiff set of crystal plasticity equations. In this case, the analytical Jacobian matrix may become unstable even due to very small difference in stress values predicted from CP and SCP models, i.e., even a very small error in stress values may lead to a very big macroscopic affine modulus M leading to an unstable analytical Jacobian matrix. Therefore, the center difference approximation is used to numerically compute the components

of the Jacobian matrix ∂ ∂ ⎯→ ⎯ E C / ˙. For instance, the component ∂ ∂ E C / 1
2 is computed by the following equation:

∂ ∂ = + - - E C E C C dh E C C dh dh ˙( ˙, ˙) ( ˙, ˙) 2 . 1 2 1 1 2 2 1 1 2 2 2 (55)
As mentioned in the beginning of the current section, we recall that the band inclination ψ and the Cauchy stress tensor in the homogeneous zone σ H are known at the end of the prescribed time increment = + t t I [ , ] n n Δ 1 . The algorithm used for the integration of the M-K governing equations over I Δ is defined by the following three main steps (Step 0 → 2). Subscript k refers to the current Newton-Raphson iteration:

• Step 0:
-Set the iteration counter k to 0. -Set C ˙1(0) and C ˙2(0) to their converged values at the previous time increment.

For ≥ k 1, do the following steps:

•

Step 1:

-Use Eq. ( 43) to compute -L k B ( 1) . -Use either the algorithm of Section 3.1 (when the conventional integration scheme is used), or the algorithm of Section 3.2 (when the DFT approach is employed), to integrate the mechanical behavior of the polycrystalline aggregate representing the band. This allows the calculation of the macroscopic Cauchy stress tensor

- σ k B ( 1)
and, in case of CP-FLD, also the macroscopic affine modulus

- M k B ( 1) . -Calculate the current iteration - f k
( 1) of the imperfection factor:

= + - - - f f ε t tL ε exp( ( ) Δ
). . Also, compute the Jacobian matrix

∂ ∂ ⎯→ ⎯ - E C C ( ˙)/ k ( 1)
, either analytically (for CP-FLD) by introducing -M k B ( 1) in Eq. ( 46), or numerically (for SCP-FLD).

-If

≤ ⎯→ ⎯ - - E C ( ˙) 10 k ( 1) L 3 2
, then the iterative procedure used to determine the jump vector C   has converged. In such a case, exit the iterative process; else, go to Step 2.

•

Step 2:

-Compute the new iteration

⎯→ ⎯ C ˙k ( ) as = + = - ⎡ ⎣ ⎢ ⎢ ∂ ∂ ⎤ ⎦ ⎥ ⎥ ⎯→ ⎯ ⎯ → ⎯ - ⎯→ ⎯ - ⎯→ ⎯ - - ⎯→ ⎯ - ⎯→ ⎯ - E E C C C C C C C ˙˙δ ˙; δ ˙( ˙) ˙• ( ˙). k k k k k k ( ) ( 1) ( 1) ( 1) ( 1) 1 ( 1) (57)
-Set ← + k k 1 and go to Step 1.

Unfortunately, the classical Newton-Raphson algorithm, defined by the previous three steps, sometimes confronts numerical difficulties. For example, when the initial guess is relatively far from the solution, the method might not converge or may converge very slowly. To avoid this shortcoming, a line search strategy leading to a global Newton-Raphson method can be applied. Further details on the implementation of this line search strategy into the Newton-Raphson iterations are provided in [START_REF] Akpama | Numerical integration of rate-independent BCC single crystal plasticity models: comparative study of two classes of numerical algorithms[END_REF].

Numerical results

It is well known that the predicted FLD is dependent on strain-rate sensitivity, plastic anisotropy, temperature and damageinduced softening, as already discussed in Section 1. In addition to these, it is also sensitive to the value of the initial imperfection factor f 0 [START_REF] Wu | On crystal plasticity FLD analysis[END_REF][START_REF] Mcginty | Application of multiscale crystal plasticity models to forming limit diagrams[END_REF]). Hence, some numerical simulations are realized using CP-FLD for deciding the value of the imperfection factor to be used for comparison with SCP-FLD results. The hardening parameters used were the same as defined in Section 3.2. The polycrystal for both homogeneous and band zones are made of 1000 single crystals possessing random initial texture. The time increment used in the different simulations of the current section is set to - 10 s 3 . Note that this time increment is much lower than that used in Section 3.2, because stress values obtained from SCP model were also used to satisfy equilibrium equation between homogeneous and band zones for FLD prediction, which requires a smaller strain step for convergence of results. Fig. 5 shows the predicted FLDs for three different values of f 0 : 0.995, 0.99 and 0.98. The results shown in Fig. 5 are consistent with the trends shown in prior literature [START_REF] Wu | On crystal plasticity FLD analysis[END_REF][START_REF] Mcginty | Application of multiscale crystal plasticity models to forming limit diagrams[END_REF].

The new spectral DFT database developed in Section 3.2 is used with the SCP model for the FLD predictions (SCP-FLD). These predictions are compared in Fig. 6 and Fig. 7 with the FLDs predicted using conventional crystal plasticity model (CP-FLD). The hardening parameters used in the different predictions of Figs. 6 and 7 are given in Table 1. For both simulations, we have fixed the initial imperfection factor f 0 to 0.99. The polycrystals used in these simulations are made of 1000 FCC single crystals possessing random initial texture. The SCP parameters used for the FLD predictions are the same as the ones defined in Section 3.2, i.e., 500 dominant DFTs were chosen for deviatoric stress components and total shear rate, and 2000 dominant DFTs for plastic spin components. As has been demonstrated in [START_REF] Knezevic | Computationally efficient database and spectral interpolation for fully plastic Taylor-type crystal plasticity calculations of face-centered cubic polycrystals[END_REF][START_REF] Knezevic | Crystal plasticity simulations using discrete Fourier transforms[END_REF], stress-strain and texture predictions improve on using a higher number of dominant DFTs. But that will also increase the computational cost. Therefore, the previous choice of number of dominant DFTs used (500 for deviatoric stress components and total shear rate, and 2000 for plastic spin components) seems to be a good compromise between accuracy and computational time in this work. If one needed to establish the FLDs for a very large number of candidate materials in designing a component, at the very least, the proposed methodology would serve to down select from the large list and perform more accurate (and expensive) computations on the shorter list. Different values of ψ Δ 0 (the step size of the initial groove orientation) were tried. We have observed that when ψ Δ 0 is equal to 1°, the SCP-FLD model allows obtaining the most optimized computations (with the minimum of CPU time and the maximum of accuracy). By analyzing Figs. 6 and7, we can easily see that the results obtained by both models (namely the SCP-FLD and CP-FLD) match very well. It can also be seen from these two figures, that the shape of FLD is different for these two sets of hardening parameters.

In Fig. 8(a), we compare the evolution of the critical strain ε 11 * , as a function of the initial band orientation ψ 0 , as predicted by the conventional (CP-FLD) and the database (SCP-FLD) approaches for two different strain-path ratios: =-ρ 0.3 and 0.3. The material parameters corresponding to Fig. 8 are the same as those used to obtain the predictions of Fig. 7. From Fig. 8(a), a similar trend can be observed for the evolution of the critical strain in the predictions from the two approaches. As already described, the necking band orientation ψ corresponding to each strain-path ratio ρ is the one minimizing the critical strain ε 11 * (among all possible initial band orientations). The plots of the necking band orientation ψ versus the strain-path ratio ρ are displayed in Fig. 8(b) for both CP-FLD and SCP-FLD computations. From Fig. 8(b), one can conclude that for positive strain-path ratios, the two approaches give exactly the same results. For this range of strain paths, the necking band orientation ψ is always equal to 0°. However, for negative strain-path ratios, the SCP model gives the limit strain at a somewhat lower necking band angle compared to the CP model. The FLDs predicted in Figs. 6 and 7 correspond to a material with random initial texture. In general, a textured material exhibits more anisotropic behavior. Fig. 9 presents a comparison of FLD predictions involving a textured material. The initial texture used in this calculation corresponds to the texture predicted by a Taylor-type model after a plane strain compression (PSC) to a true strain of -0.5. The hardening parameters of the material are taken the same as in Fig. 7. The results from conventional CP-FLD computations reported in Fig. 9(a) reveal that, in the negative range of strain-path ratios, the shape and the level of the FLD predicted for the textured material are almost the same as those with random initial texture. However, the two FLDs differ significantly in the positive range of strain-path ratios, especially in the vicinity of equibiaxial stretching. This observation is consistent with other reports in the literature [START_REF] Barlat | Crystallographic texture, anisotropic yield surfaces and forming limits of sheet metals[END_REF][START_REF] Wu | Effect of cube texture on sheet metal formability[END_REF][START_REF] Yoshida | The effects of texture on formability of aluminum alloy sheets[END_REF][START_REF] Signorelli | Predictions of forming limit diagrams using a rate-dependent polycrystal self-consistent plasticity model[END_REF]. It has been found that the initial and induced plastic anisotropy mainly influence the yield surface shape and its evolution during plastic deformation. The yield surface shape has, in turn, a tremendous effect on the forming limit diagrams. For the biaxial stretching range, it has been demonstrated that a decrease in the sharpness of the yield locus in the equibiaxial stretching is a factor promoting larger limit strains [START_REF] Neale | Limit strain predictions for strain-rate sensitive anisotropic sheets[END_REF]. The good agreement between the CP-FLD and SCP-FLD in Fig. 9(b) provides support for the claims made earlier that the databases developed here are indeed independent of the initial texture.

Since the accuracy of both CP-FLD and SCP-FLD models is nearly the same, we compare the CPU times required by these techniques to obtain the FLDs. We remind the reader that this is the main motivation behind this work. We have observed in Section 3.2 that, compared to the CP model, the SCP model is two orders of magnitude faster with same degree of accuracy. The accuracy of the SCP model can be increased by using a larger number of dominant DFTs. It should be noted that even using a small number of dominant DFTs we are able to achieve similar accuracy as of the CP model. Similar observation w.r.t. accuracy between CP-FLD and SCP-FLD is observed in Figs. 6,7 and 9(b). Various parallelization techniques can be implemented in the conventional crystal plasticity code for FLD predictions. These might include solving single crystal equations in parallel for polycrystal or solving for different strain-path ratios (ρ) in parallel. Although these would reduce the required computational time, the requirement for computational resources increases. In present study, parallelization was implemented only for CP-FLD simulations in solving the single crystal equations. CP-FLD simulations were done on a cluster at Georgia Tech using 4 processors of clock speed 2.5 GHz each and 12 GB memory. The total CPU time (sum of time taken by each processor) for predicting FLD using CP model on an average was ∼500 h. As can be seen, parallelization was essential in this case to get the solutions in a reasonable timeframe. However, the SCP-FLD simulations were performed on the same cluster using a single processor of clock speed 2.5 GHz and 12 GB memory. The total CPU time for FLD prediction using SCP model on an average is ∼20 h. Thus, the computational time (CPU time) required for predicting FLD reduces dramatically when using SCP model, resulting in about 96% savings. Indeed, one can even predict the FLD using the proposed SCP-FLD approach on a normal desktop computer/laptop. This is a significant advantage compared to the high computational resources requirement for CP-FLD. On a normal desktop computer/laptop with Intel Core i3 processor of clock speed 3.6 GHz and 4 GB memory, the SCP-FLD simulation takes ∼12 h. The gain in computational speed between CP and SCP model in Section 3.2 was observed as two orders of magnitude, which is slightly lower for FLD prediction using these two models. This lower than expected gain in the computational efficiency is due to the additional calculations involved in the FLD prediction, such as those requiring the solution of the equilibrium equation between homogeneous and band zones using Newton-Raphson iteration and subroutine calls for calculating average polycrystal stress response using crystal plasticity thousands of times. Nevertheless, the SCP model is still able to predict FLD by using only 4% of the computational time required for predicting FLD by conventional crystal plasticity model.

The studies published in literature are generally not focused on reducing the computational time for FLD prediction; hence, it is difficult to find previous reports that mention the computational time for prediction of full FLD along with reporting all the details of material and simulation parameters to be able to make a fair comparison. Therefore, due to the absence of all relevant details, it is not possible to compare our approach in terms of computational time to other numerically efficient algorithms for predicting FLD, like the one discussed in [START_REF] Schwindt | Experimental and Visco-Plastic Self-consistent evaluation of forming limit diagrams for anisotropic sheet metals: an efficient and robust implementation of the M-K model[END_REF] where a robust numerical procedure based on direct MK-VPSC approach is proposed. Spectral database approach is aimed at reduced-order modeling of crystal plasticity theories for general applications (within some idealizations and approximations as discussed) and essentially retains all features of a crystal plasticity theory, while being fast in evaluations. Most other reduced-order models are narrower in their scope, i.e., specific to an application. Given the highly nonmonotonic strain paths simulated in arriving at FLDs, the present case study presents the toughest validation of the DFT approach to date and proposes a computationally efficient SCP-FLD tool.

Conclusions

In this paper, we have developed a highly computationally efficient tool for predicting Forming Limit Diagrams (FLDs) based on rate-dependent crystal plasticity models using the DFT-based spectral database. Significant improvements in the accuracy of the crystal plasticity computations using DFT approach (SCP) were obtained by using new spectral DFT database, which implemented improvements like discretization of uniform grid of orientation and deformation mode space to one degree, separate databases for all six components of the deviatoric stress and use of trigonometric interpolation for recovering values of function not on grid points. It was also demonstrated that the SCP model is able to predict the FLD by using only 4% of the computational time compared to conventional crystal plasticity model leading to 96% savings in computational cost. Moreover, the SCP-FLD tool allows the prediction of full FLD, implementing polycrystalline anisotropy on a normal desktop/laptop containing fast processor (e.g., Intel Core i3) in about half a day time. Thus this computationally efficient tool leads to major savings of computational cost in FLD predictions. Interestingly, it has been demonstrated by [START_REF] Knezevic | A high-performance computational framework for fast crystal plasticity simulations[END_REF][START_REF] Mihaila | Three orders of magnitude improved efficiency with high-performance spectral crystal plasticity on GPU platforms[END_REF][START_REF] Savage | Computer implementations of iterative and non-iterative crystal plasticity solvers on high performance graphics hardware[END_REF] that there can be a further gain in the computational speed of SCP compared to CP, if the computations are performed on a specialized computer hardware that utilizes graphics-processing units (GPU). Thus, there is clear potential to further speed up the SCP-FLD tool presented here. (1 1 1) (1 1 1) (1 1 1) (1 1 1) (1 1 1)

α 7 8 9 1 0 1 1 1 2 ⎯ → ⎯⎯⎯ m 2 α 0 [1 01] [0 1 1] [1 1 0] [0 1 1] [1 1 0] [1 0 1] ⎯→ ⎯ n 3 α 0 (1 1 1)
(1 1 1) (1 1 1) (1 1 1) (1 1 1) (1 1 1)

Tensors P α 0 and A 0 α involved in Eqs. ( 29) and ( 36

  that define respectively the normal to the slip plane and the slip direction for the slip system α in the current configuration. Further details on vectors are given in Appendix A.

  (1) and (2) allows us to obtain the following relation between the strain rate tensor D (the symmetric part of L) and tensor S:

Fig. 1 .

 1 Fig. 1. M-K analysis for a metal sheet: (a) Initial geometry of the sheet; (b) Current geometry of the sheet.

  frame of D in which all crystal plasticity calculations using database can be performed and converted back to sample reference frame when needed. Hence, g e ( ,{ }) j p

Fig. 2 .

 2 Fig. 2. Variation of components D D , 1 2 and D 3 as a function of θ.

Fig. 3 .

 3 Fig. 3. Comparison of predictions from SCP model using new (one_deg_grid) and old (three_deg_grid) database against corresponding predictions from conventional CP Taylor-type model for uniaxial tension of OFHC Copper: (a) pole figures of deformed texture, and (b) different components (S11 and S 33 ) of deviatoric stress vs. strain curves.

Fig. 4 .

 4 Fig. 4. Comparison of predictions from SCP model using new (one_deg_grid) and old (three_deg_grid) database against corresponding predictions from conventional CP Taylor-type model for equibiaxial tension of OFHC Copper: (a) pole figures of deformed texture, and (b) different components (S 11 and S 33 ) of deviatoric stress vs. strain curves.

B

  can be extracted from the inverse of the macroscopic affine modulus -M 1 (see Eq. (

Fig. 5 .

 5 Fig. 5. CP-FLD predictions for three different values of initial imperfection factor.

Fig. 6 .

 6 Fig. 6. Predicted FLDs using conventional crystal plasticity model (CP-FLD) and spectral crystal plasticity model using DFT database (SCP-FLD) for the following hardening parameters: = h 180 MPa 0

Fig. 7 .

 7 Fig. 7. Predicted FLDs using conventional crystal plasticity model (CP-FLD) and spectral crystal plasticity model using DFT database (SCP-FLD) for the following hardening parameters: = h 500 MPa 0

Fig. 9 .

 9 Fig. 9. Impact of the initial texture: (a) Comparison of the CP-FLDs predicted with an initially random texture and a textured material (texture corresponding to plane strain compression to a true strain of -0.5 on the initially random texture using a Taylor-type model); (b) Comparison of the CP-FLD and SCP-FLD predictions for the textured material.

Fig. 8 .

 8 Fig. 8. (a) Evolution of the critical strain ε 11* versus the initial band orientation ψ 0 as predicted by both CP-FLD and SCP-FLD models for two particular strain-path ratios ( =-ρ 0.3 and 0.3); (b) Evolution of the necking band orientation ψ as a function of the strain-path ratio ρ as predicted by both CP-FLD and SCP-FLD models.

  be constant in the isoclinic configuration and are enumerated in Table A.1.Table A.1The numbering of the slip systems of an FCC single crystal.

  ), respectively, are the symmetric and skew-symmetric parts of the Schmid tensor

  The set of constitutive equations for the single crystal presented in Section 2.1 is integrated, over a typical time increment =

	3. Numerical aspects
	3.1. Conventional integration of the crystal plasticity constitutive equations
	t t I [ , n n Δ	+	] 1 (where	t n	+	1	-	t n
									B
			1	2		2 1			B r	11 B 22 12 B	1 11 H 1 12 H	H 2 12 H 2 22	(28)

Table 1

 1 Hardening parameters used in the simulations of Figs. 6 and 7.

		h 0 (MPa)	s s (MPa)	a	s 0 (MPa)	m
	Fig. 6	180	148	2.25	16	0.01
	Fig. 7	500	230	2.8	50	0.01

Acknowledgements

Akash Gupta and Surya R. Kalidindi gratefully acknowledge support received for this work from the Office of Naval Research under the award number N00014-15-1-2478 (Dr. William M. Mullins, program manager). Akash Gupta also acknowledges support received from management of Tata Consultancy Services for this work.

Appendix A. Slip systems for FCC single crystals

The slip direction ⎯ → ⎯⎯⎯