domain: spi phys

text:
Vikram Phalke
email: vikramphalke040@gmail.com

Samuel Forest
email: samuel.forest@minesparis.psl.eu

Hyung-Jun Chang
email: hyung-jun.chang@safrangroup.com

Tonya Rose
email: tonya.rose@safrangroup.com

Arjen Roos
email: arjen.roos@safrangroup.com

Dark energy in crystals. Prediction of stored energy in polycrystalline aggregates

Keywords: Stored energy, Crystal plasticity, Taylor-Quinney coefficient, Single crystals, Polycrystals, Strain gradient plasticity, GND, Micromorphic crystal plasticity

During the plastic deformation of metallic materials, part of expended mechanical energy diffuses as heat. The fraction of plastic work converted into heat is called the Taylor-Quinney Coefficient (TQC), which is often assumed to be a constant parameter of about 0.9. The remaining portion of the plastic work is called stored energy. The stored energy is known as the main driving force for dynamic or static recovery and recrystallization. Therefore, numerical predictions and experimental measurements of the stored energy and TQC are essential to optimize thermomechanical material processing. An adequate prediction of the stored energy and the TQC using existing crystal plasticity models in line with the experimental measurements remains a challenging problem. In this work, a thermodynamic class of crystal plasticity models is used to predict the stored energy and TQC of copper and aluminum single crystals. Then, the numerical stored energy predictions are extended to polycrystalline austenitic steel 316L and compared with the experimental measurements from the literature. An adhoc factor is introduced in the numerical expression of stored energy in order to compensate for the difference with the experimental measurement. To this end, the contributions of statistically stored dislocations (SSDs) and geometrically necessary dislocations (GNDs) for the stored energy prediction are analyzed to understand the physical origin of the ad-hoc factor. The contribution of GNDs to stored energy and enhanced hardening is accounted for by means of a strain gradient plasticity model. The present systematic finite element crystal plasticity simulations also include specific interface conditions at grain boundaries. The presented computational analysis indicates that, compared to the experiment, there remains dark energy in the evaluation of energy storage as predicted by the proposed thermodynamically consistent crystal plasticity framework.

Introduction

During the plastic deformation of metallic materials, the fraction of plastic work converted into heat is called the Taylor-Quinney Coefficient (TQC), which is often assumed to be a constant parameter of about 0.9. The remaining portion of the plastic work is called stored energy. The stored energy is the main driving force for dynamic or static recovery and recrystallization. The prediction of the stored energy and TQC is important to understand the plastic deformation and subsequent recovery and recrystallization mechanisms [START_REF] Rollett | Recrystallization and Related Annealing Phenomena[END_REF].

Different techniques used in experimental measurements of stored energy are summarized in [START_REF] Bever | The stored energy of cold work[END_REF]. Extensive experimental work has been performed in the past to measure the stored energy in metallic materials, for instance in [START_REF] Bever | The stored energy of cold work[END_REF][START_REF] Bailey | The dislocation density, flow stress and stored energy in deformed polycrystalline copper[END_REF][START_REF] Williams | The stored energy of copper deformed at 24°C[END_REF][START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF][START_REF] Oliferuk | Effect of the grain size on the rate of energy storage during the tensile deformation of an austenitic steel[END_REF]. Ravichandran et al. [START_REF] Ravichandran | On the conversion of plastic work into heat during high strain rate deformation[END_REF] investigated the TQC evolution in aluminum 2024-T3 alloy and α-titanium. They found that the TQC is a function of strain but not of strain rate in the case of aluminum 2024-T3 alloy. On the other hand, TQC is strongly dependent on strain rate in α-titanium. Rittel et al. [START_REF] Rittel | On the dynamically stored energy of cold work in pure single crystal and polycrystalline copper[END_REF] investigated the thermo-mechanical response of single and polycrystalline copper at low and high strain rates. They observed that the TQC linearly increases with strain rate and remains considerably lower than the classic value of 0.9. Rittel et al. [START_REF] Rittel | The dependence of the Taylor-Quinney coefficient on the dynamic loading mode[END_REF] recently measured the TQCs for 7 different metals and alloys, namely, Ti6Al4V, commercially pure titanium, Al5086, Al2024, 304L, 1020 steel, and maraging C300 under tension, compression, and dominant shear loading. A significant deviation of TQC from 0.9 was observed. Except for commercially pure titanium, identical TQCs were measured under tension, compression, and shear loading. In commercially pure titanium, the presence of twinning activity in compression and shear loading results in different TQC for different deformation modes. Global measurements of plasticity induced heating are usually performed although intergranular temperature field measurements are becoming possible, at least under quasi-static loading conditions [START_REF] Charkaluk | Microplasticity in Polycrystals: A Thermomechanical Experimental Perspective[END_REF].

Several attempts have been made to investigate the effect of grain size on stored energy, see for instance [START_REF] Williams | The stored energy of copper deformed at 24°C[END_REF][START_REF] Oliferuk | Effect of the grain size on the rate of energy storage during the tensile deformation of an austenitic steel[END_REF][START_REF] Baker | The effect of grain size on the stored energy of cold work as a function of strain for polycrystalline nickel[END_REF] which show a decrease in stored energy for increasing grain size. Williams [START_REF] Williams | The stored energy of copper deformed at 24°C[END_REF] found that an increase in grain size by a factor of ten decreases the stored energy by 8% at a strain of 0.3. Recently, Rittel et al. [START_REF] Rittel | On the dynamically stored energy of cold work in pure single crystal and polycrystalline copper[END_REF] measured the stored energy of single crystal and polycrystalline copper in the strain rate range of 1000 -8000 s -1 undergoing compressive loading. Higher stored energy was found in polycrystals than single crystals.

An adequate prediction of stored energy and TQC in line with the experimental measurements using numerical simulations remains a challenging issue. A hierarchy of material models for crystalline solids is presented in [START_REF] Mesarovic | Mesoscale Models. From Micro-Physics to Macro-Interpretation[END_REF] where discrete (molecular and dislocation dynamics) methods are used to inform continuum crystal plasticity models. Analytic expressions used to predict the stored energy may not consider all the mechanisms responsible for stored energy. Nietofuentes et al. [START_REF] Nieto-Fuentes | On a dislocation-based constitutive model and dynamic thermomechanical considerations[END_REF] introduced an adhoc factor in the analytic expression of stored energy to represent the experimentally measured stored energy adequately. This requires calibration of this factor for several dynamically loaded materials. A well-known analytic expression for the stored energy taken as a function of dislocation density, can be found in [START_REF] Bever | The stored energy of cold work[END_REF][START_REF] Bailey | The dislocation density, flow stress and stored energy in deformed polycrystalline copper[END_REF]. Bailey [START_REF] Bailey | The dislocation density, flow stress and stored energy in deformed polycrystalline copper[END_REF] measured the stored energy of polycrystalline silver and copper under tensile loading. The analytic expression used by Bailey [START_REF] Bailey | The dislocation density, flow stress and stored energy in deformed polycrystalline copper[END_REF] for the prediction of stored energy volume density is given by E s = (τc) 2 µ , where the critical resolved shear stress τ c = 0.5µb √ ρ avg depends on the average dislocation density ρ avg . Furthermore, he plotted the experimentally measured stored energy, E s , as a function of (τc) 2 µ . The results for both metals lie close to the straight line of equation E s ≃ 7.7 (τc) 2 µ ≃ 2µb 2 ρ avg . Crystal plasticity modeling is a convenient tool to determine the stored energy and evolution of the material microstructure in thermo-mechanical processes. Hakansson et al. [START_REF] Håkansson | Prediction of stored energy in polycrystalline materials during cyclic loading[END_REF] investigated the effect of initial texture on the stored energy in polycrystalline austenitic steel 316L using a rate-dependent crystal plasticity model for large deformations formulated within a thermodynamic framework. A finite strain formulation of crystal plasticity in a thermodynamically consistent manner is used by Jafari et al. [START_REF] Jafari | A finite-deformation dislocation density-based crystal viscoplasticity constitutive model for calculating the stored deformation energy[END_REF] to predict the stored energy in single and polycrystalline aluminum under tensile loading and bi-crystal aluminum under compressive loading. Note that the reliability of the crystal plasticity model is crucial to predict the stored energy. Models of plastic deformation fall into two categories. On the one hand, phenomenological models are frequently used in crystal plasticity modeling. However, these models do not give a physical connection with the microscopic mechanisms of plastic deformation. Physics-based crystal plasticity models can overcome this limitation by introducing physically motivated internal variables such as dislocation density in the constitutive framework. The constitutive equations must be derived in a thermodynamically consistent manner [START_REF] Cailletaud | An overview of the anatomy of crystal plasticity models[END_REF]. In particular, it is necessary to check the positivity of the internal entropy production rate. Failure to do so can violate the first and second laws of thermodynamics [START_REF] Ottosen | Thermodynamic framework for constitutive modeling[END_REF]. Thermodynamic consistency of physics-based models has received less attention compared to phenomenological models [START_REF] Wu | Thermodynamic considerations on a class of dislocation-based constitutive models[END_REF].

Numerically, dislocation density-based models are often used to describe material hardening. The obtained total dislocation density is further used to calculate stored energy. The energy stored in copper single crystals under tensile loading using discrete dislocation plasticity was prdicted by Benzerga et al. [START_REF] Benzerga | The stored energy of cold work: Predictions from discrete dislocation plasticity[END_REF]. They showed that the dislocation distribution influences the stored energy along with the total dislocation density. Moreover, they found TQCs in the range of 0.75-0.95 up to the strain of 10% and strain rate of 100 s -1 . A reliable dislocation density-based model is necessary to adequately represent the total dislocation density and consequently the stored energy according to [START_REF] Nieto-Fuentes | On a dislocation-based constitutive model and dynamic thermomechanical considerations[END_REF]. Kositski et al. [START_REF] Kositski | Employing molecular dynamics to shed light on the microstructural origins of the Taylor-Quinney coefficient[END_REF] used molecular dynamics simulations to show that the dislocation accumulation mechanism is not the only mechanism responsible for stored energy. They include grain boundary evolution as an additional mechanism responsible for energy storage. The difficulty in providing a thermodynamic framework for plasticity stems from the fact that evolution of dislocation structures is often viewed as analogous to turbulence [START_REF] Kubin | Dislocations, Mesoscale Simulations and Plastic Flow. Oxford Series on Materials Modelling[END_REF]. This is reflected by the fact that the complex Kocks-Mecking-Teodosiu dislocation density evolution equations are not easily incorporated in a complete thermodynamic setting involving free energy density potential and a dissipation function [START_REF] Maugin | Internal variables and dissipative structures[END_REF][START_REF] Maugin | The Thermomechanics of Nonlinear Irreversible Behaviors[END_REF]. One reason is the presence of phenomenological interaction and hardening matrices with complex evolution [START_REF] Franciosi | Crystal hardening and the issue of uniqueness[END_REF].

According to Ashby's terminology, dislocation densities can be split into densities of statistically stored dislocations (SSD) and densities of geometrically necessary dislocations (GND) [START_REF] Ashby | The deformation of plastically non-homogeneous alloys[END_REF]. The latter concept is closely related to the dislocation density tensor or Nye's tensor which has been included in a thermodynamically consistent constitutive framework in [START_REF] Steinmann | Views on multiplicative elastoplasticity and the continuum theory of dislocations[END_REF][START_REF] Shizawa | A thermodynamical theory of gradient elastoplasticity with dislocation density tensor. I: Fundamentals[END_REF]. GND densities can be estimated from measurements of crystal lattice curvature by means of EBSD methods [START_REF] Pai | Study of grain boundary orientation gradients through combined experiments and strain gradient crystal plasticity modeling[END_REF]. Most of the numerical work on the prediction of stored energy considers the contribution of SSDs only. Studies on the prediction of stored energy considering the contribution of GNDs along with SSDs are relatively rare. Proper evaluation of GND densities is possible by means of strain gradient plasticity (SGP) models incorporating the dislocation density tensor or the lattice curvature tensor [START_REF] Mesarovic | Mesoscale Models. From Micro-Physics to Macro-Interpretation[END_REF]. Estimates for stored energy values from crystal plasticity analyses is currently used for the prediction of annealing and recrystallization phenomena [START_REF] Chatterjee1 | Evaluating the influence of deformation variables on dynamic recrystallization behavior using a crystal plasticity model[END_REF] but also crack initiation and propagation under fatigue loading conditions, see [START_REF] Xu | Microstructural fracture mechanics: Stored energy density at fatigue cracks[END_REF] for HCP alloys and [START_REF] Shen | Stored energy density research on the fatigue crack initiation at twin boundary and life prediction of Inconel718 superalloy[END_REF] for Inconel 718. Confrontations of such predictions to experimental measurements available in the literature for several FCC metals and alloys have been the subject of little analysis to date.

The aim of the present work is to evaluate the ability of a class of thermodynamic formulations of standard and strain gradient crystal plasticity models to predict the stored energy and TQC for some single and poly-crystalline FCC metallic materials. At first, the necessity to ensure the positivity of the entropy production rate is emphasized. In addition, it is shown that, for some physics-based models, this can be ensured by imposing constraints on some model parameters. Then the stored energy and TQC are predicted using the standard crystal plasticity model for copper and aluminum single crystals. Next, the application is made to polycrystalline aggregates. As a first application to polycrystalline simulations, the effect of mesh size and grain morphology on the stress-strain response and energy evalution is studied briefly. The stored energy predicted for polycrystalline austenitic steel 316L is compared with the experimentally measured stored energy from the literature. A constant ad-hoc factor, ξ, in the expression of stored energy, akin to the one introduced in [START_REF] Nieto-Fuentes | On a dislocation-based constitutive model and dynamic thermomechanical considerations[END_REF], is treated as a fitting parameter so that the predicted stored energy adequately represents the experimental measurements. The next application deals with energy stored in annealed Inconel 718 deformed at high strain rates (2000 s -1 ) and a larger temperature range (293 K-1073 K). Very high temperatures are excluded in order to avoid recrystallization phenomena, not addressed in this work. In addition, the predicted stored energy considering contribution of both SSDs and GNDs is compared to that obtained by considering the contribution of SSDs only. Finally, the contribution of GNDs along with SSDs is discussed in the prediction of stored energy using two different grain boundary conditions, i.e. intermediate and microhard.

The novel aspects addressed in the present work are the following: A consistent thermodynamic framework accommodating Kocks-Mecking-Teodosiu evolution equations and estimates of GND densities, the evaluation of the stored energy density for this class of models for various single and polycrystalline FCC metals and alloys, and comparison of the results with several measurements from the literature. It will be shown that, in the absence of the corrective ad-hoc factor, the models systematically underestimate the amount of stored energy compared to available measurements from the literature, hence the name dark energy adopted in the title. The study is limited to quasi-static processes, including rate-dependent constitutive laws, even though some high strain rates will be imposed at some stage for comparison with literature results. The materials retained in the analysis are chosen either because crystal plasticity models or measurements of stored energy are available.

The paper is organized as follows. The constitutive setting is described in Section 2. Section 3 is devoted to analyzing the positivity of the dissipation rate to ensure thermodynamic consistency. The summary of the constitutive equations used for the numerical prediction of stored energy and TQC is given in section 4. In section 5, single crystals aluminum and copper simulations are performed to predict the stored energy and evolution of TQCs. In addition, this section provides the prediction of stored energy and evolution of TQC in two polycrystalline FCC engineering metallic materials, with appropriate discussion of the results. Conclusions follow in section 6.

The following notations are employed in this contribution: Underlined A and under-waved bold A ∼ characters are used to denote first-order and second-order tensors, respectively. The transpose, inverse and time derivative are represented as

A ∼ T , A ∼ -1
and Ȧ ∼ . Simple and double contractions are understood in the sense a • b = a i b j and A ∼ : B ∼ = A ij B ij . Gradient and divergence operators are defined with respect to Lagrange coordinates, Grad , Div , and with respect to Euler coordinates, grad , div .

Constitutive framework for thermo-elasto-viscoplastic single crystals

Kinematics

The material points are defined by their position vector X in the reference configuration Ω 0 and their current position vector x in the configuration Ω t . The large deformation framework of thermo-elasto-viscoplasticity is based on the multiplicative decomposition of the total deformation gradient F ∼ into a recoverable thermo-elastic part F ∼ the and a plastic part F ∼ p following [START_REF] Bertram | Finite thermoplasticity based on isomorphisms[END_REF][START_REF] Ristinmaa | Thermodynamic format and heat generation of isotropic hardening plasticity[END_REF][START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF]:

F ∼ = F ∼ the • F ∼ p . ( 1
)
The spatial velocity gradient and the plastic deformation rate are defined as

ℓ ∼ = Ḟ ∼ • F ∼ -1 , L ∼ p = Ḟ ∼ p • F ∼ p-1 . ( 2
)
The volume mass densities with respect to the reference configuration, the intermediate configuration, and the current configuration are called ρ 0 , ρ # and ρ, respectively, and given by

J = det(F ∼ ) = ρ 0 ρ , J the = det(F ∼ the ) = ρ # ρ , J p = det(F ∼ p ) = ρ 0 ρ # . ( 3
)
In dense metallic crystals considered in this work, plastic deformation is isochoric such that

J p = det F ∼ p = 1, J the = detF ∼ the = J = det F ∼ , ρ # = ρ 0 . ( 4
) However, J p can be different from one in the case of compressible plasticity. This situation was studied for ductile fracture of porous single crystals in [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]. Moreover, the thermo-elastic Green-Lagrange strain tensor E ∼ the is introduced as follows:

E ∼ the = 1 2 [(F ∼ the ) T • (F ∼ the ) -1 ∼ ], (5)
with 1 ∼ denoting the second order identity tensor.

Definition of stress tensors

The second Piola-Kirchhoff stress tensor Π ∼ e with respect to the intermediate configuration is defined as

Π ∼ e = J the (F ∼ the ) -1 • σ ∼ • (F ∼ the ) -T . ( 6
)
where σ ∼ is the Cauchy stress tensor. The tensor Π ∼ e is related to the thermo-elastic strain tensor E ∼ the by

Π ∼ e = Λ ≈ : (E ∼ the -α(T -T 0 )1 ∼ ). (7)
where Λ ≈ is the fourth-order tensor of elastic moduli and α is the thermal expansion coefficient for cubic materials. The first Piola-Kirchhoff stress tensor Π ∼ p related to the Cauchy stress tensor σ ∼ is given by

Π ∼ p = Jσ ∼ F ∼ -T . ( 8
)
Moreover, the Mandel stress tensor Π ∼ M which is work-conjugate to L ∼ p in the power of internal forces, can be defined with respect to the intermediate configuration as follows:

Π ∼ M = J the (F ∼ the ) T • σ ∼ • (F ∼ the ) -T . ( 9
)

Flow rule

A rate-dependent overstress-type flow rule is adopted, based on a Schmid yield function defined as

f r = |τ r | -τ r c , ( 10
) involving the resolved shear stress τ r on the slip system r, which is the driving force to trigger plastic slip, and the corresponding critical resolved shear stress τ r c . The resolved shear stress τ r on slip system r is defined as

τ r = Π ∼ M : (m r ⊗ n r ), (11)
where m r is the slip direction and n r is the slip plane normal vector in the lattice frame.

The slip rate γr on each slip system r is then given by the following rate-dependent flow rule

γr = f r K m sign(τ r ), (12)
with Macauley brackets < • > denoting the positive part of •, and K and m are the viscosity parameters. The high values of power m and low values of K lead to an almost rate-independent elasto-plastic behavior in a given strain rate range. Plastic deformation rate is the result of slip processes on N distinct slip systems, according to the flow rule:

L ∼ p = N r=1 γr (m r ⊗ n r ). ( 13
)
The cumulative plastic strain variable γ cum is introduced as

γ cum = t 0 N r=1 | γr |dt. (14)

Dislocation density-based hardening model

Phenomenological flow rules are frequently used in crystal plasticity modeling, following the original works [START_REF] Asaro | Crystal plasticity[END_REF]. They share the drawback that the material state is only described in terms of critical resolved shear stress and not in terms of lattice defect population such as dislocation densities. This limitation can be overcome by resorting to physics-based crystal plasticity models. The latter provide a strong physical relationship with the microscopic mechanisms of plastic deformation by introducing internal variables such as dislocation density in the constitutive framework. The dislocation density-based models have better predictability compared to the phenomenological models over a wide range of strain, strain rates, and temperatures [START_REF] Čebron | Stored energy predictions from dislocation-based hardening models and hardness measurements for tensile-deformed commercial purity copper[END_REF].

In the present work strain hardening behavior relies on a dislocation density-based hardening model, which takes dislocation interactions into account. Following [START_REF] Kubin | Modeling dislocation storage rates and mean free paths in face-centered cubic crystals[END_REF], the critical resolved shear stress τ r c is directly related to the scalar dislocation densities as follows:

τ r c = τ 0 + µ N u=1 h ru ϱ u , ( 15
)
where τ 0 is the lattice friction stress, h ru is the interaction matrix describing interaction between dislocations, ϱ u = ρ u b 2 is the non-dimensional dislocation density for slip system u, and ρ u is the usual dislocation density, i.e. the length of dislocation lines per unit volume with b being the norm of the dislocation Burgers vector. The following equation gives the evolution of the dislocation density [START_REF] Teodosiu | Finite element simulation of the large elastoplastic deformation of multi-crystals[END_REF]:

ρr = | γr |   N u=1 b ru ϱ u κ c -d c ϱ r   . ( 16
)
called the Kocks-Mecking-Teodosiu evolution equation. The dislocation interaction is described by the matrix b ru , κ c is a constant material parameter proportional to the number of obstacles crossed by a dislocation before being immobilized, and y c = bd c is the critical distance controlling the annihilation of dislocations with opposite signs. The structure of the matrices h ru and b ru can be found in [START_REF] Ling | A reduced micromorphic single crystal plasticity model at finite deformations. application to strain localization and void growth in ductile metals[END_REF]. The total non-dimensional dislocation density is defined as follows:

ϱ = N r=1 ϱ r 0 + t 0 N r=1 ρr dt, ( 17
)
where ϱ r 0 is the initial non-dimensional dislocation density for slip system r.

Micromorphic crystal plasticity model

The previous standard crystal plasticity must now be enhanced to evaluate the effect of GNDs on hardening and energy storage. For that purpose, the reduced order micromorphic crystal plasticity model proposed by [START_REF] Ling | A reduced micromorphic single crystal plasticity model at finite deformations. application to strain localization and void growth in ductile metals[END_REF] is used as an alternative to strain gradient plasticity models based on the full dislocation density tensor, due to its computational efficiency. The predictions of this model have been compared with strain gradient plasticity formulations involving the dislocation density tensor in [START_REF] Phalke | Modeling size effects in microwire torsion: A comparison between a lagrange multiplier-based and a CurlF p gradient crystal plasticity model[END_REF]. Results under monotonic tensile loading and torsion show that both models have similar predictions regarding size-dependent hardening due to GNDs, at a lower computational price for the reduced micromorphic model. The micromorphic approach is an efficient method for the numerical implementation of strain gradient crystal plasticity as discussed in [START_REF] Sansour | A formulation for the micromorphic continuum at finite inelastic strains[END_REF][START_REF] Cordero | Size effects in generalised continuum crystal plasticity for two-phase laminates[END_REF]. The main equations of the latter model are summarized in this section following [START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF].

According to the micromorphic approach, the material points possess two types of degrees of freedom: The displacement vector u (X , t) = x -X and the micromorphic scalar microslip variable γ χ (X , t). The gradients of the degrees of freedom with respect to the reference configuration are called

∂u ∂X = Grad u , K = ∂γ χ ∂X = Grad γ χ . ( 18
)
The static balance equations and Neumann boundary conditions expressed with respect to the reference configuration are as follows:

Div Π ∼ p = 0 and Div M -S = 0, ∀ X ⊂ Ω 0 , (19)
T = Π ∼ p • N and M = M • N , ∀ X ⊂ ∂Ω 0 , ( 20
)
where S and M are generalized stresses conjugate to microslip and its gradient, respectively, M is the generalized surface traction and N the outward unit normal vector at the boundary of the reference body. Note that the Cauchy stress tensor is still symmetric in the proposed micromorphic theory.

The microslip variable γ χ is related to the cumulative plastic strain γ cum , see Eq. ( 14), via the relative plastic slip variable e p defined as

e p := γ cum -γ χ . ( 21
)

General thermodynamic framework

The material under consideration is characterized by the coupled thermo-mechanical Helmholtz free energy density function defined in terms of the thermo-elastic strain tensor E ∼ the , the relative plastic slip e p , the gradient of the microslip variable K , temperature T and the set of internal hardening variables ζ as follows:

Ψ = Ψ(E ∼ the , T, ζ, e p , K ). ( 22
)
The Helmholtz free energy density function is assumed to take the form:

ρ 0 Ψ(E ∼ the , T, ζ, e p , K ) = 1 2 E ∼ the : Λ ≈ : E ∼ the + 1 2 H χ e 2 p + 1 2 K • A ∼ • K + ρ 0 C ε (T -T 0 ) -T log T T 0 + (T -T 0 )P ∼ : E ∼ the + ρ 0 Ψ ζ (T, ζ) ( 23
)
where Ψ ζ denotes the part of the free energy density associated with the internal variables and will be detailed in Section 3. The expression involves the tensor Λ ≈ of the elastic moduli and generalized micromorphic moduli H χ and A ∼ . The higher order moduli A ∼ have the physical dimension Pa.m 2 ≡N which means that they contain constitutive internal lengths that can be obtained for instance after dividing them by the shear modulus. The modulus H χ is used as a numerical penalty parameter ensuring that e p almost vanishes, i.e. that the micromorphic variable is equal to the cumulative plastic slip [START_REF] Scherer | Lagrange multiplier based vs micromorphic gradient-enhanced rate-(in)dependent crystal plasticity modelling and simulation[END_REF]. The tensor P ∼ is related to thermal expansion and T 0 is a chosen reference temperature. Quadratic potentials are therefore adopted for the elastic/reversible contributions. The logarithmic temperature contribution is discussed for instance in [START_REF] Bertram | On the introduction of thermoplasticity[END_REF]. It is such that the heat capacity under constant strain conditions, defined as C ε = -T ∂ 2 Ψ/∂ 2 T is constant in the temperature range [START_REF] Ristinmaa | Thermodynamic format and heat generation of isotropic hardening plasticity[END_REF]. The linearization of this expression with respect to T -T 0 leads to the usual quadratic contribution in the infinitesimal case.

The second thermodynamic principle, or entropy imbalance, is written in the form of the local Clausius-Duhem inequality following [START_REF] Germain | Continuum Thermodynamics[END_REF]:

Jσ ∼ : ℓ ∼ -ρ 0 ( Ψ + η Ṫ ) -Q • Grad T T ≥ 0 ( 24
)
where η is the mass density of entropy and Q is the heat conduction vector with respect to the reference configuration, given by Q = J the F ∼ -1 • q and Grad T the Lagrangian gradient of temperature. Expanding the time derivative of the free energy density function and the power of internal forces leads to the following form of the Clausius-Duhem inequality

Π ∼ e -ρ 0 ∂Ψ ∂E ∼ the : Ė ∼ the -S + ρ 0 ∂Ψ ∂e p ėp + M -ρ 0 ∂Ψ ∂K • K + S γcum + Π ∼ M : L ∼ p -ρ 0 η + ∂Ψ ∂T Ṫ -ρ 0 ∂Ψ ∂ζ ζ -Q • Grad T T ≥ 0. (25)
The three first contributions related to standard and micromorphic elasticity are assumed to be reversible which leads to the following state laws:

Π ∼ e = ρ 0 ∂Ψ ∂E ∼ the , S = -ρ 0 ∂Ψ ∂e p , M = ρ 0 ∂Ψ ∂K ( 26
)
The usual relation between temperature and entropy is implemented:

η = - ∂Ψ ∂T ( 27
)
The thermodynamic forces associated with the internal variables are defined as

X = ρ 0 ∂Ψ ∂ζ . ( 28
)
As a result, the dissipation rate, which restricts the material flow and hardening rules in connection with the yield condition, is given by

D := Π ∼ M : L ∼ p + S γcum -X ζ -Q • Grad T T ≥ 0. ( 29
)
The thermodynamic forces associated with the arguments of the Helmholtz free energy function are derived from the potential [START_REF] Maugin | The Thermomechanics of Nonlinear Irreversible Behaviors[END_REF]:

Π ∼ e = Λ ≈ : (E ∼ the -E ∼ th ), S = -H χ e p = -H χ (γ cum -γ χ ), M = A ∼ • K . ( 30
)
In the micromorphic theory, the introduced yield function ( 10) is enhanced in the following way

f r = |τ r | -(τ r c -S) = |τ r | -(τ r c -DivM ).
(31) after considering the second micromorphic balance law [START_REF] Benzerga | The stored energy of cold work: Predictions from discrete dislocation plasticity[END_REF]. For cubic materials, the second-order tensor takes the simple form A ∼ = A1 ∼ , where A is the generalized modulus assumed to be uniform. The viscoplastic flow rule ( 12) is still valid for each slip system in the micromorphic case except that the yield function ( 31) is used instead of (10).

Thermodynamic consistency of the crystal plasticity constitutive equations

In this section, several existing hardening models are evaluated from the perspective of thermodynamic consistency. Following the general thermodynamic framework settled in subsection 2.6, the analysis concentrates on the contribution of the set ζ of internal variables to stored energy Ψ ζ (T, ζ) and to the dissipation rate. The micromorphic part of the model is dropped in this section for simplicity.

The non-negative dissipation rate D from Eq. ( 24) can be split into two parts: Non-negative mechanical dissipation rate D m ≥ 0 and non-negative heat condition term

D th ≥ 0 such that D = D m + D th ≥ 0 with D m = Jσ ∼ : ℓ ∼ -ρ 0 ( Ψ + η Ṫ ) and D th = -Q • Grad T T .
Positivity of D m and D th are required independently. The thermal inequality D th ≥ 0 is called Fourier inequality which states that heat must flow from hot to cold regions. Finally, the dissipation rate reads

D m = Π ∼ M : L ∼ p -X ζ ≥ 0, ( 32
)
where the Mandel stress tensor Π ∼ M is defined by Eq. ( 9). In crystal plasticity involving N slip systems, at least N internal variables ζ r must be introduced. Their meaning will depend on the specific material model. The corresponding part of the free energy potential in Eq. ( 23) is additively decomposed as

Ψ ζ (T, ζ r ) = N r=1 Ψ r ζ (T, ζ r ), and X r = ρ 0 ∂Ψ r (ζ r ) ∂ζ r , (33)
so that the dissipation rate in Eq. ( 32) can be further expressed as

D m = N r=1 τ r γr - N r=1 X r ζr , with Π ∼ M : L ∼ p = N r=1 τ r γr , ( 34
)
where X r is the thermodynamic force associated with each internal hardening variable ζ r . Any acceptable set of constitutive relations must fulfill the dissipation inequality. The amount of dissipated energy is fundamentally dependent on the choice of internal variables ζ r , the expression of the free energy potential and the evolution equations for internal variables. Various approaches ensuring fulfillment of the dissipation inequality can be found in [START_REF] Ottosen | Thermodynamic framework for constitutive modeling[END_REF][START_REF] Lemaitre | Mechanics of Solid Materials[END_REF][START_REF] Besson | Non-Linear Mechanics of Materials[END_REF]. In the first approach called direct approach, some evolution laws are assumed for ζ r . The positivity of the dissipation rate is then checked at each time step of the computation. Positivity may then be ensured only for special values of material parameters and specific loading paths. A second approach is called Onsager's approach in which linear relationships between the rate of internal hardening variables ζr and thermodynamic forces X r are assumed [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF][START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF]. The linear relationship cannot be used for general plasticity and visco-plasticity as realistic materials show nonlinearity. The third approach is called potential approach and generalizes the Onsager method to nonlinear theory. In that case, evolution equations for internal variables and plastic flow are derived from a dissipation potential function [START_REF] Maugin | The Thermomechanics of Nonlinear Irreversible Behaviors[END_REF][START_REF] Germain | Continuum Thermodynamics[END_REF]. If this function is convex with respect to its argument, positivity of dissipation rate is ensured for a given range of material parameters and for all loading paths. The existence of the internal energy potential is a requirement of the first thermodynamic principle. The situation is different for the dissipation potential whose existence is not ensured by the laws of nature. For complex sets of evolution equations, such a dissipation potential may not exist, as it will become apparent in the next examples.

The existence of a dissipation potential is not necessary, but it can ease the development of thermodynamically consistent models, i.e. models that automatically fulfill the second law of thermodynamics. These potentials are not always provided in the literature, but they are useful to finally evaluate the TQC.

Example of phenomenological crystal plasticity model

At first, a thermodynamically consistent phenomenological model is considered, inspired by [START_REF] Besson | Non-Linear Mechanics of Materials[END_REF][START_REF] Busso | On the selection of active slip systems in crystal plasticity[END_REF], in which the contribution of the internal hardening variable ζ r to the free energy is assumed to have the quadratic form:

ρ 0 Ψ r (ζ r ) = 1 2 Qζ r N s=1 h rs ζ s , ( 35
)
where Q is a material parameter and h rs an interaction matrix. The thermodynamic force X r associated with the internal hardening variable ζ r therefore is

X r = ρ 0 ∂Ψ r (ζ r ) ∂ζ r = Q N s=1 h rs ζ s . ( 36
)
In this model, the existence of a convex dissipation potential Ω(Π ∼ M , X r ) is assumed from which the flow rule and internal variable evolution equations are derived as:

Ḟ ∼ p • F ∼ p-1 = ∂Ω ∂Π ∼ M , ζr = - ∂Ω ∂X r . ( 37
)
In addition, the dissipation potential is assumed to be of the power law form:

Ω(Π ∼ M , τ r c ) = K m + 1 N r=1 f r (Π ∼ M , τ r c ) K m+1 , ( 38
)
where τ r c is the critical resolved shear stress and f r is the Schmid-type yield function [START_REF] Charkaluk | Microplasticity in Polycrystals: A Thermomechanical Experimental Perspective[END_REF]. Based on the dissipation potential in Eq. [START_REF] Kubin | Modeling dislocation storage rates and mean free paths in face-centered cubic crystals[END_REF], the flow and hardening rules are derived

Ḟ ∼ p • F ∼ p-1 = ∂Ω ∂Π ∼ M = N r=1 ∂Ω ∂f r ∂f r ∂Π ∼ M = N r=1 γr (m r ⊗ n r ), ζr = - ∂Ω ∂τ r c = vr , ( 39
)
where v r is the cumulative slip variable for slip system r, and γr is the slip rate on slip system r computed as follows:

γr = ∂Ω ∂f r sign(τ r ) = vr sign(τ r ), vr = | γr |. ( 40
)
The dissipation rate in Eq. ( 34) can be written as:

D m = N r=1 τ r γr -Q N r=1 N s=1 h rs v s | γr |, with ζr = | γr | = N r=1 |τ r | -τ 0 + Q N s=1 h rs v s + τ 0 | γr |. (41)
Taking ζ r = v r , the second term in the brackets is exactly the value of the critical resolved shear stress. With the yield function [START_REF] Charkaluk | Microplasticity in Polycrystals: A Thermomechanical Experimental Perspective[END_REF] and for τ 0 ≥ 0, it appears that each term is either zero for negative values of the yield function, or positive otherwise. This proves that the dissipation inequality is always fulfilled for this model.

The constitutive formulation presented above is rather simple as ζ r = v r . It is more realistic to use a saturating variable ζ r instead of the cumulative slip v r , which is not a fully satisfactory internal variable. In phenomenological models, ζ r is then reminiscent of a dislocation density-like variable. Materials generally display nonlinear hardening behavior which motivates a nonlinear isotropic hardening rule for ζ r such as

ζ r = 1 -exp(-Bv r ), with vr = | γr |, (42)
where B ≥ 0 is a material parameter. The contribution of ζ r to the free energy is assumed to be of the same form as in Eq. [START_REF] Ling | An elastoviscoplastic model for porous single crystals at finite strains and its assessment based on unit cell simulations[END_REF]. Then the thermodynamic force associated with the internal hardening variable becomes

X r = ρ 0 ∂Ψ r (ζ r ) ∂ζ r = Q N s=1 h rs (1 -exp(-Bv s )). ( 43
)
In addition, the nonlinear evolution of the isotropic hardening variable τ r c is taken of the form

τ r c = τ 0 + X r = τ 0 + Q N s=1 h rs (1 -exp(-Bv s )). ( 44
)
The dissipation rate in [START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF] then reads

D m = N r=1 |τ r | -Q N s=1 h rs (1 -exp(-Bv s ))(B exp(-Bv s )) | γr |, (45)
which can be further expressed as follows:

D m = N r=1 |τ r | -τ 0 + Q N s=1 h rs (1 -exp(-Bv s )) + τ 0 + X r (1 -B exp(-Bv r )) | γr |. (46)
Again, in a way similar to the analysis of ( 41), it appears that each term in the sum [START_REF] Germain | Continuum Thermodynamics[END_REF] is either zero or positive. This is due to the fact the first two terms in the brackets are nothing but the yield function f r itself and the two last contributions are positive for suitable values of material parameters.

This class of phenomenological single crystal models has been applied to engineering applications, see for instance [START_REF] Coudon | A stochastic approach applied to directionally solidified turbine blades[END_REF]. It has the advantage that no check of positive dissipation at each material point and each time step is required.

Physics-based crystal plasticity models

Phenomenological models have the advantage that free energy functions are explicitly postulated; sometimes, a dissipation potential is also proposed, as done in the previous section. But the models may be too simple compared to experimental results. In contrast, physics-based internal variables related to deformation mechanisms may display better prediction capabilities. The dislocation densities ρ r on the slip systems r are considered as physically relevant internal variables, for instance in [START_REF] Kubin | Modeling dislocation storage rates and mean free paths in face-centered cubic crystals[END_REF][START_REF] Knezevic | Multiscale modeling of microstructure-property relationships of polycrystalline metals during thermo-mechanical deformation[END_REF][START_REF] Bronkhorst | Structural representation of additively manufactured 316l austenitic stainless steel[END_REF]. But explicit expressions of the free energy function are not provided in the literature in most cases. In contrast, Lieou et al. [START_REF] Lieou | Thermodynamic theory of crystal plasticity: Formulation and application to polycrystal fcc copper[END_REF] derived the evolution of dislocation density from energetic and entropic considerations alone with constraints of the first and second laws of thermodynamics. The formulation is based on the Langer-Bouchbinder-Lookman thermodynamic dislocation theory proposed in [START_REF] Langer | Thermodynamic theory of dislocation-mediated plasticity[END_REF][START_REF] Langer | Statistical thermodynamics of strain hardening in polycrystalline solids[END_REF]. Based on the energetic considerations alone, the evolution of dislocation density is proportional to plastic work rate and not only to the plastic slip rate as in conventional theories. These model formulations however differ from the continuum thermodynamics framework adopted in the present work which is inherited from [START_REF] Maugin | The Thermomechanics of Nonlinear Irreversible Behaviors[END_REF][START_REF] Germain | Continuum Thermodynamics[END_REF][START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]. In this section, it is shown how the latter thermodynamic framework can accommodate a class of physics-based crystal plasticity models used in this paper.

The hardening laws introduced in subsection 2.4 are analyzed in the present section from the perspective of thermodynamic consistency. Following the work of [START_REF] Abrivard | Phase field modelling of grain boundary motion driven by curvature and stored energy gradients. part I: theory and numerical implementation[END_REF], the contribution of the internal hardening variable ζ r to the free energy function is assumed to have the quadratic form:

ρ 0 Ψ r (ζ r ) = 1 2 µξ(ζ r ) 2 , with ζ r = N s=1 h rs ϱ s , ( 47
)
where ξ is an ad-hoc factor or material parameter which is expected to be of the order of unity. Then, the thermodynamic force X r associated with ζ r is given by

X r = ρ 0 ∂Ψ r (ζ r ) ∂ζ r = µξζ r . ( 48
)
Moreover, the evolution of the critical resolved shear stress, regarded as the isotropic hardening variable, τ r c , is given by Eq. ( 15). This means that the thermodynamic force X r = ξ(τ r c -τ 0 ) in the present formulation. The evolution of non-dimensional dislocation density ϱ r (= ρ r b 2 ) is described by Eq. ( 16). After inserting the expression Eq. ( 16) for ρr , the expression for dissipation related to the internal hardening variable becomes

N r=1 X r ζr = 1 2 µξ N r=1 N s=1 h rs ρr = 1 2 µξ N r=1 N s=1 h rs N u=1 b ru ϱ u κ c -d c ϱ r | γr | . ( 49
)
Inserting previous equation in [START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF] gives the dissipation rate of the form:

D m = N r=1 |τ r | - 1 2 µξ N s=1 h rs N u=1 b ru ϱ u κ c -d c ϱ r | γr |. (50)
Enforcing the positivity of the dissipation rate in (50) is not possible in all situations. It must therefore be checked at each material point and each time step in numerical simulations for given values of the material parameters.

Let us consider two special cases for which positivity of the dissipation rate can be ensured under any loading conditions by imposing constraints on the model parameters.

Case 1: According to Taylor's assumption of isotropic interaction between slip systems, all components of the dislocation interaction matrix h rs and b rs , which determine the Taylor stress τ r c and dislocation mean free path, respectively, are unity such that

D m = N r=1 |τ r | - 1 2 µξ N u=1 ϱ u κ c -d c ϱ r | γr |, (51)
The dissipation rate in the previous equation can also be expressed by using τ r c as follows:

D m = N r=1 |τ r | -τ r c + τ 0 + µ N u=1 ϱ u - µξ 2κ c N u=1 ϱ u + 1 2 µξd c ϱ r | γr |, (52)
which further gives

D m = N r=1 (|τ r | -τ r c ) + τ 0 + µ N u=1 ϱ u 1 - ξ 2κ c + 1 2 µξd c ϱ r | γr |. ( 53
)
Each individual term of the sum in the previous equation is zero or positive respectively when | γr | is 0 or the corresponding term in the brackets is positive. The first term in the brackets is the Schmid yield function Eq. ( 10) which is positive when | γr | ̸ = 0. The third term can be made positive by imposing constraints on the model parameter κ c such that κ c ≥ ξ/2. Besides, τ 0 ≥ 0 and the last term of the previous equation is also positive. As a result, under these circumstances, the dissipation rate is positive for all loading conditions of the material points and the thermodynamic consistency is automatically fulfilled. The thermodynamic consistency of the Kocks-Mecking-Teodosiu evolution equation is rarely discussed in the literature. In the case of isotropic latent hardening, we have derived a constraint on the dimensionless parameters κ c and ξ from the second principle. The parameter ξ is an "ignorance" parameter for the evaluation of stored energy from the dislocation densities while κ c is the "ignorance" parameter for the multiplication rate of dislocations. Their default order of magnitude is 1, which does not violate the found constraint κ c ≥ ξ/2. The parameter κ c much be significantly greater than zero to avoid too high multiplication rate, which again fulfills the constraint. These parameters however remain phenomenological.

Case 2: It is assumed that the dislocation interaction matrices h rs and b rs are diagonal matrices with all diagonal components equal to 1. This simplified case, different from Taylor isotropic interaction, turns out to be relevant in the case of nickel-based single crystal superalloys [START_REF] Méric | Single crystal modeling for structural calculations. Part 1: Model presentation[END_REF]. It follows that the evolution of the isotropic hardening variable τ r c becomes

τ r c = τ 0 + µ √ ϱ r , ( 54
) and ( 50) can be expressed as follows:

D m = N r=1 (|τ r | -τ r c ) + τ 0 + µ √ ϱ r 1 - ξ 2κ c + 1 2 µξd c ϱ r | γr |. (55)
The same reasoning as in the previous case shows that when κ c ≥ ξ/2, positivity of the dissipation rate is satisfied for all loading conditions of the material point.

In the present work, the general case of dissipation rate presented in Eq. ( 50) is considered. The dissipation rate will be evaluated at each material point and at each time increment in order first to check its positivity and, second, to evaluate the dissipated energy by time integration.

Expressions for stored energy, temperature rise and Taylor-Quinney coefficient for numerical simulations

In this section, the formulas for estimating the stored energy, Taylor-Quinney coefficient, and temperature evolution in the context of standard and micromorphic crystal plasticity models are summarized.

Expressions for standard crystal plasticity models

Stored energy density

The expression for the stored energy based on the proposed thermodynamically consistent framework, cf. Eq. ( 47), for the standard crystal plasticity model is given by

E c = ρ 0 N r=1 Ψ r (ζ r ) = 1 2 µξ N r=1 (ζ r ) 2 , with ζ r = N s=1 h rs ϱ s . ( 56
)
The volume averaged stored energy over the volume V of any considered material body is

E c avg = 1 V V E c dV. ( 57
)

Total non-dimensional dislocation density

The total non-dimensional dislocation density is defined by Eq. ( 17). The volume averaged non-dimensional dislocation density over a material volume V is computed as

ϱ avg = 1 V V ϱ dV. ( 58
)

Temperature evolution under adiabatic conditions

Expressions for the temperature evolution under adiabatic conditions for standard and micromorphic crystal plasticity are derived in [START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF]. The temperature evolution for the standard crystal plasticity model is given by

Ṫ = Π ∼ M : Ḟ ∼ p F ∼ p-1 -X ζ ρC ε . ( 59
)
Inserting ( 34) and ( 49) into the previous equation gives

Ṫ = N r=1 τ r γr -1 2 µξ N r=1 N s=1 h rs ρr ρC ε . ( 60
)
Note that this expression does not contain the thermo-elastic contribution since it is negligible compared to temperature rise induced by plasticity [START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF][START_REF] Russo | Thermomechanics of cosserat medium: Modeling adiabatic shear bands in metals[END_REF].

Taylor-Quinney coefficient (TQC)

The integral form of TQC is taken as

β int = t 0 N r=1 τ r γr dt -1 2 µξ t 0 N r=1 N s=1 h rs ρr dt t 0 N r=1 τ r γr dt . ( 61
)
The volume averaged Taylor-Quinney coefficient over a material volume V is computed as

β avg = 1 V V β int dV. ( 62
)

Expressions in the context of micromorphic crystal plasticity model

Stored energy density

Combining the general expression [START_REF] Maugin | The Thermomechanics of Nonlinear Irreversible Behaviors[END_REF] of the micromorphic energy density with the dislocation density based contribution ( 33) and ( 47), the stored energy for the micromorphic crystal plasticity model is taken as

E s = 1 2 µξ N r=1 (ζ r ) 2 + 1 2 H χ e 2 p + 1 2 AK • K . ( 63
)
The two last terms represent estimates of the stored energy due to GNDs. The norm of the microslip gradient vector is regarded as a fair approximation of the norm of the dislocation density tensor although these quantities are not identical, as discussed in [START_REF] Phalke | Modeling size effects in microwire torsion: A comparison between a lagrange multiplier-based and a CurlF p gradient crystal plasticity model[END_REF] and also in [START_REF] Abatour | Toward robust scalar-based gradient plasticity modeling and simulation at finite deformations[END_REF] regarding the limit of this energy contribution when H χ → ∞, i.e. when the micromorphic model converges toward the strain gradient plasticity model, strictly speaking.

The volume averaged stored energy over any material volume V is given by

E s avg = 1 V V E s dV. ( 64
)

Temperature evolution under adiabatic conditions

The temperature evolution for the micromorphic crystal plasticity model is given by

Ṫ = N r=1 τ r γr + S γcum -1 2 µξ N r=1 N s=1 h rs ρr ρC ε . ( 65
)
It contains an additional contribution S γcum in the dissipation rate compared to the standard expression [START_REF] Russo | Thermomechanics of cosserat medium: Modeling adiabatic shear bands in metals[END_REF]. Again, the thermoelastic contribution is dropped so that the adiabatic temperature can only increase during deformation. This is used as a check of the positivity of the dissipation rate which coincides with the numerator in [START_REF] Abrivard | A coupled crystal plasticity -phase field formulation to describe microstructural evolution in polycrystalline aggregates during recrystallisation[END_REF].

Taylor-Quinney coefficient (TQC)

The integral form of TQC is given by

β int = t 0 N r=1 (τ r γr + S γcum ) dt -1 2 µξ t 0 N r=1 N s=1 h rs ρr dt t 0 N r=1 (τ r γr + S γcum ) dt . ( 66
)
The volume averaged Taylor-Quinney coefficient over any material volume is also given by Eq. ( 62).

Crystal plasticity Finite element simulation results

Single crystals simulations

In a preliminary study, simulations are performed for aluminum and copper single crystals subjected to tensile loading to predict the evolution of stored energy and TQC. Finite element (FE) simulations are performed based on the constitutive framework of standard crystal plasticity model presented in section 2. A simplified geometry of 0.06 mm × 0.06 mm × 0.06 mm cube is used, as in [START_REF] Alankar | A dislocation-density-based 3D crystal plasticity model for pure aluminum[END_REF], having 8 C3D8 type elements which are 8 node linear brick elements. The applied boundary conditions are presented in Fig. 1. They lead to homogeneous straining of single crystals with high symmetry with respect to the cube. Deviation from high symmetry orientations leads to heterogeneous deformation.

<001> and <111> crystal orientations

Two high symmetry orientations of the single crystal considered are <001> and <111> such that the [001] or [111] crystal direction align with the loading direction, respectively. In FCC crystals, the crystallographic slip can occur on the 12 {111}<110> slip systems. The experimental stress-strain responses from the work [START_REF] Hosford | Tensile deformation of aluminum single crystals at low temperatures[END_REF] for aluminum single crystals and [START_REF] Takeuchi | Work hardening of copper single crystals with multiple glide orientations[END_REF] for copper single crystals are used to calibrate material parameters. A misorientation of < 1 • is applied from the tensile axis as in the experimental tests of [START_REF] Hosford | Tensile deformation of aluminum single crystals at low temperatures[END_REF] and [START_REF] Takeuchi | Work hardening of copper single crystals with multiple glide orientations[END_REF], as considered also in [START_REF] Abrivard | A coupled crystal plasticity -phase field formulation to describe microstructural evolution in polycrystalline aggregates during recrystallisation[END_REF]. All their experiments were performed at room temperature and applied strain rates of 7.5 × 10 -5 s -1 and 10 -3 s -1 for aluminum and copper single crystals, respectively. The material parameters given in Table 1 are essentially taken from the literature but the viscosity and Kocks-Mecking parameters were further calibrated for better Table 1: Material constants and fitted material parameters for the simulation of single crystalline aluminum under tensile loading. The interaction matrix coefficients (h 0 -h 5 ) are taken from [START_REF] Kubin | Modeling dislocation storage rates and mean free paths in face-centered cubic crystals[END_REF].

C11 C12 C44 ρ Eq. (60)

Cε

Eq. ( 60) τ0 (fitted) Eq. ( 15) b 108 GPa 61.3 GPa 28.5 GPa 2700 kg m -3 900 Jkg -1 K -1 0.6 MPa 0.286 nm µ ϱ r 0 (fitted) Eq. ( 17)

K (fitted)
Eq. ( 12) m (fitted) Eq. ( 12) κc (fitted) Eq. ( 16) dc (fitted) Eq. ( 16) 27 GPa 8 × 10 -10 0.5 MPa.s

Cε

Eq. ( 60) τ0 (fitted) Eq. ( 15) b 170 GPa 124 GPa 75 GPa 8960 kgm -3 385 Jkg -1 K -1 0.5 MPa 0.257 nm µ ϱ r 0 (fitted) Eq. ( 17)

K (fitted)
Eq. ( 12) m (fitted) Eq. ( 12) κc (fitted) Eq. ( 16) dc (fitted) Eq. ( 16) 41 GPa 8 × 10 -10 0.5 MPa.s description of the curves at various strain rates. In particular the coefficients of the hardening matrix h rs are taken from [START_REF] Kubin | Modeling dislocation storage rates and mean free paths in face-centered cubic crystals[END_REF], with parameter notations from [START_REF] Gérard | Modeling of latent hardening produced by complex loading paths in fcc alloys[END_REF]. The dislocation density-based hardening model captures the stress-strain responses in reasonable agreement with the experimentally measured responses for <001> and <111> crystal orientations as shown in Fig. 2. The material hardening parameter κ c mainly governs the initial slope of the stress-strain curve, while parameter d c controls the saturation. Numerical values of the material constants and fitted material parameters for the single crystals aluminum and copper are summarized in Table 1 and 2, respectively. The material parameters are fitted using <001> crystal orientation (Fig. 2). The tensile axes <001> and <111> are oriented for multi-slip with 8 and 6 equally favored slip systems, respectively. The initial hardening rate of the <111> crystal orientation is higher than the <001> crystal orientation. The initial dislocation density ρ r 0 (= ϱ r 0 /b 2 ) is assumed to be 10 10 m -2 for both single crystals and chosen the same for all slip systems. The total dislocation density evolution for single crystals aluminum and copper is presented in Fig. 3a and 3b, respectively. These figures show that the dislocation density increases rapidly and saturates with further increase in plastic strain. The stored energy is predicted using a thermodynamically consistent formulation of the standard crystal plasticity model, cf. Eq. [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]. Fig. 4 shows that the stored energy is strongly orientation dependent. Similar to the dislocation density evolution, the stored energy increases rapidly and saturates with further increase in plastic strain. The variation of TQCs for aluminum and copper single crystals is shown in Fig. 5a and 5b, respectively. The TQCs predicted for the aluminum and copper single crystals remain around 0.95 which is above the commonly measured value of 0.9.

Note that only the increase in dislocation density and stored energy is plotted on the curves of Fig. 3a, 3b and all other figures in the sequel. The initial value of total dislocation density is at least 1000 times smaller than the values reached during tensile straining according to the model after 0.05 strain. Accordingly, the curves are essentially unaffected by adding initial values.

<123> crystal orientation

Rittel et al. [START_REF] Rittel | On the dynamically stored energy of cold work in pure single crystal and polycrystalline copper[END_REF] measured the TQC for single crystal copper of <123> orientation under compressive loading subjected to quasi-static and high strain rates deformation. Fig. 4: Predicted stored energy E s using the thermodynamically consistent formulation (Eq. ( 56)) at a temperature of 298 K for (a) aluminum and (b) copper single crystals (<001> and <111> orientations). Taylor-Quinney coefficent

< 111 > < 001 > (b)
Fig. 5: Predicted Taylor-Quinney coefficient (β int ) using the thermodynamic formulation (Eq. ( 61)) for the single crystals (<001> and <111>) at a temperature of 298 K for (a) aluminum and (b) copper.

The average TQC values were in the range of rather low values 0.2 -0.3 in quasistatic deformation in contrast to 0.65 -0.85 for high-strain rate deformation tests. To compare the experimentally measured TQC values, simulations are performed with <123> crystal orientation subjected to compressive loading with direction [123] parallel to the loading direction. The material parameters used in the simulations are given in Table 2, the same as previously. The experimentally measured stress-strain responses and corresponding predicted responses are shown in Fig. 6a. The viscosity parameters used capture the strain rate sensitivity observed in the experimental work. Note that temperature evolves under adiabatic conditions but the impact on the material parameters is not taken into account here, in constrast to our recent work on adiabatic shear banding [START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF]. A discrepancy can be observed in the overall stress-strain responses after 0.08 strain, which may be due to thermal softening not accounted for here. The <123> crystal orientation triggers a single slip at the initial deformation stage, but multiple slip systems are activated in a later stage.

The predicted dislocation density evolution is presented in Fig. 6b. The corresponding stored energy and TQCs at two strain rates (0.1 s -1 and 3000 s -1 ) are shown in Fig. 7a and7b, respectively. The high strain rate deformation test leads to more stored energy than low strain rate tests (Fig. 7a). It is found that the predicted TQCs are higher than the values measured by Rittel et al. [START_REF] Rittel | On the dynamically stored energy of cold work in pure single crystal and polycrystalline copper[END_REF] (Fig. 7b). One possible reason for this discrepancy is that the total dislocation density predicted by the model is Taylor-Quinney coefficent

3000 s -1 0.1 s -1 (b)
Fig. 7: Evolution of the predicted (a) stored energy according to Eq. ( 56) and (b) Taylor-Quinney coefficient (Eq. ( 61)) for single crystal copper (<123>) with ξ = 1 under compressive loading using boundary conditions of Fig. 1.

significantly lower than the actual experimental values but these were not measured. The total dislocation density necessary to predict the stored energy and consequently the TQC needs to be at least of the order of 10 15 m -2 . In contrast, the dislocation density obtained by numerical simulations is of the order of 10 14 m -2 according to Fig. 6b. However, higher dislocation densities will lead to higher and unrealistic stress values, according to formula [START_REF] Jafari | A finite-deformation dislocation density-based crystal viscoplasticity constitutive model for calculating the stored deformation energy[END_REF]. Another possible explanation is that energy formula (47)) is not precise enough and more sophisticated functions of dislocation densities are required. This expression of the stored energy itself may not take into account all existing mechanisms responsible for the stored energy, as suggested in [START_REF] Williams | The stored energy of copper deformed at 24°C[END_REF].

Simulation of polycrystalline aggregates

In this section, polycrystalline simulations are performed to predict the stored energy and TQC in a few metals and alloys. First of all, a study of the effect of mesh size and grain morphology on the volume averaged stress-strain responses over the whole FE model and local fields is performed. Next, the prediction of stored energy and TQC is carried out using the standard crystal plasticity model. Finally, the impact of GND densities by means of the micromorphic crystal plasticity model on the predicted stored energy is evaluated.

Effect of mesh size and grain morphology

In general, representative volume elements (RVE) with cubic meshes are used to predict the texture evolution and the global stress-strain responses, for instance in [START_REF] Alankar | A dislocation-density-based 3D crystal plasticity model for pure aluminum[END_REF][START_REF] Kalidindi | Crystallographic texture evolution in bulk deformation processing of fcc metals[END_REF]. These authors found that the global stress-strain curves using RVEs having cubic meshes are well represented. However, the local fields inside grains were not considered. More realistic polycrystalline microstructures can be represented using tetrakaidecahedra [START_REF] Mika | Effects of grain interaction on deformation in polycrystals[END_REF][START_REF] Mika | Polycrystal plasticity modeling of intracrystalline boundary textures[END_REF], or Voronoï tessellation as in [START_REF] Barbe | Intergranular and intragranular behavior of polycrystalline aggregates. part 1: F.e. model[END_REF][START_REF] Diard | Evaluation of finite element based analysis of 3d multicrystalline aggregates plasticity: Application to crystal plasticity model identification and the study of stress and strain fields near grain boundaries[END_REF]. In some studies, experimentally determined 3D microstructures are incorporated in crystal plasticity modeling so that the local fields and global stress-strain responses can be compared to experimental results, see for example [START_REF] Musienko | Three-dimensional finite element simulation of a polycrystalline copper specimen[END_REF][START_REF] Proudhon | Coupling diffraction contrast tomography with the finite element method[END_REF].

RVEs with cubic meshes

A brief study of the effect of mesh size and grain morphology on the stress-strain behavior is performed using the standard crystal plasticity model according to section 2. At first, three RVEs of 0.3 mm × 0.3 mm × 0.3 mm, discretized with structured mesh using C3D20R elements, which are 20 node quadratic brick elements with reduced integration, are considered in Fig. 8. Each RVE consists of 64 grains assigned with random orientations for isotropic crystallographic texture. To study the mesh size each grain of the RVE is discretized with 1 element (Fig. 8a), 8 elements (Fig. 8b), and 27 elements (Fig. 8c). The material parameters used in this polycrystalline study are given for copper in Table 2. The applied boundary conditions to the RVE are presented in Fig. 1, corresponding to mixed homogeneous boundary conditions.

The predicted volume averaged stress-strain responses using three RVEs are displayed in Fig. 9a. The stress-strain curve obtained using the RVE with 27 elements per grain displays a significantly softer response compared to RVEs with 1 and 8 elements per grain (keeping 64 grains per RVE). This confirms the need for sufficient discretization of intragranular fields even for the prediction of overall curves.

In addition, tests are performed to study the effect of the number of grains on the stress-strain behavior using RVEs having 64, 512, and 1000 grains (keeping 1 elements per grain). Fig. 9b shows that there is no significant effect of number of grains on the stress-strain behavior for these loading conditions and isotropic texture.

RVEs generated by Voronoï tessellations

In addition to cuboidal grains, RVEs of 0.3 mm × 0.3 mm × 0.3 mm are generated by Voronoï tessellations, see Fig. 10. To study the effect of mesh size, an RVE containing 64 grains is meshed with two different mesh sizes, namely the coarse (1774 nodes, 8274 elements) and fine (16402 nodes, 85251 elements) using C3D4R elements, which are 4-node linear tetrahedral elements with reduced integration. The cumulative plastic strain fields are shown in Fig. 12. The fine meshed RVE captures the heterogeneity of the local fields in contrast to almost homogeneous field predicted by coarse meshed RVE. This suggests that fine enough mesh size within each grain is required to predict the heterogeneity of the local fields. Moreover, the predicted stress-curves are shown in Fig. 13a. The coarsely meshed RVE displays stiffer response compared to the finely meshed RVE, keeping the same total number of grains.

In addition to the effect of mesh size, the effect of grain number on the averaged stress-strain response is studied. The three RVEs considered for the study contain 64, 125, and 200 grains displayed in Fig. 11. The predicted averaged stress-strain curves of Fig. 13b show that hundred grains are sufficient for a stabilized overall response for the considered tensile test.

Prediction of stored energy in engineering polycrystals

Two engineering alloys are considered for the prediction of stored energy in polycrystals. Firstly, stored energy and TQC are computed for polycrystalline austenitic steel 316L under tensile loading at room temperature at low strain rate using the standard crystal plasticity model of section 2. This material has been chosen because the predicted stored energy can be compared with the experimental data taken from [START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF].

Secondly, the method is applied to Inconel 718 under compression loading at high strain rates and high temperatures because material parameters are available with their temperature dependence from [START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF].

Stored energy and TQC in polycrystalline austenitic steel 316L

The material parameters τ 0 , ϱ r 0 , κ c , d c , K, and m are calibrated against the experimental stress-strain data of [START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF] for an average grain size of 80 µm. The coefficients of the interaction matrix (h 0 -h 5 ) are adopted from [START_REF] Hure | Intergranular stress distributions in polycrystalline aggregates of irradiated stainless steel[END_REF]. The initial dislocation density ρ r 0 (= ϱ r 0 /b 2 ) is assumed to be 10 10 m -2 and is the same for all slip systems. The initial dislocation density is an essential information. In combination with the friction stress it sets the initial yield stress and also the initial hardening rate via the Kocks-Mecking-Teodosiu equation. It follows that it is not a free parameter that we could calibrate with impunity to better estimate the amount of stored energy. The RVE of 0.27 mm × 0.27 mm × 0.27 mm is used to describe the polycrystalline austenitic steel. It is generated using Voronoï tessellation having 64 grains assigned with random orientations. The applied boundary conditions to the RVE are those of Fig. 1. The fitted averaged stress-strain response against the experimental response is shown in Fig. 14a. The material constants and fitted material parameters are summarized in Table 3.

Estimation of stored energy and TQC is performed using the expressions ( 56) and ( 61), respectively. The predicted volume average stored energy (Eq. ( 57)) predicted by the standard crystal plasticity model using ξ = 1 is plotted against the experimental data of [START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF] in Fig. 14b. The figure clearly shows that the predicted amount of stored energy is significantly underestimated by a factor of about 10. Similar observation has been made for single crystal copper (<123> crystal orientation) simulations where numerically predicted stored energy and, consequently, the TQC were underestimated compared to experimental results, see section 5.1. Again, this raises questions about the actual value of dislocation density in the samples (not measured) and on the reliability of the simple energy evaluation [START_REF] Langer | Thermodynamic theory of dislocation-mediated plasticity[END_REF]. In addition, the temperature rise measured in the experiments can be compared with the numerical predictions. Oliferuk et al. [START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF] report a temperature evolution from 25.5°C to 28.5°C during tensile tests on austenitic steel at a strain rate of 0.002 per second, up to 7% prescribed strain. The Table 3: Material parameters used for the numerical simulation of polycrystalline austenitic steel 316L under tensile loading at room temperature. The elastic constants and dislocation interaction coefficients are taken from [START_REF] Hure | Intergranular stress distributions in polycrystalline aggregates of irradiated stainless steel[END_REF].

C11 C12 C44 ρ Eq. (60)

Cε

Eq. ( 60) τ0 (fitted) Eq. ( 15) 199 GPa 136 GPa 105 GPa 7965 kgm

-3 532 Jkg -1 K -1 80 MPa b µ ϱ r 0 (fitted)
Eq. ( 17)

K (fitted)
Eq. ( 12) m (fitted) Eq. ( 12) κc (fitted) Eq. ( 16 simulation predicts a temperature rise of 15 • C in the same conditions which is 5 times higher. Two reasons can be advocated to explain the discrepancy. First, adiabatic conditions set an upper bound for the temperature rise. These conditions are questionable for the considered quasi-static tests. Second, the underestimation of stored energy by the model leads to higher local temperatures. This pleads for more systematic and accurate measurements combined with constitutive model development and thermomechanical structural simulations following early attempts such as [START_REF] Chrysochoos | Plastic and dissipated work and stored energy[END_REF].

The ad-hoc parameter ξ can be increased so as to improve the agreement between computational and experimental results. The predicted volume average stored energy with ξ = 10 is shown in Fig. 14b. With such a high value of parameter ξ, the predicted stored energy is in line with the experimental measurements. Moreover, the evolution of TQC with strain is shown for ξ = 1, 10 in Fig. 15. The predicted TQC with ξ = 1 varies between 0.95 to 0.98. The values of TQC measured by Oliferuk et al. [START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF] were in the range 0.58-0.75. With ξ = 10, the predicted TQC is in the range of 0.55 to 0.70. The expression of stored energy given by Bailey [START_REF] Bailey | The dislocation density, flow stress and stored energy in deformed polycrystalline copper[END_REF] for both polycrystalline silver and copper gives an ad-hoc factor of the order of 2 as explained in section 1. In contrast, in this work, for polycrystalline austenitic steel based on the numerical simulations of polycrystals, we find an ad-hoc factor of the order of 10.

Stored energy and TQC in Inconel 718

After emphasizing the importance of an ad-hoc factor in the prediction of stored energy and TQC, the application is made to polycrystalline Inconel 718 undergoing higher strain rates and higher temperature compression loading. The predicted stress-strain 56), ( 57)) and comparison against the experimental measurements from the work of [START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF]. responses and corresponding experimental responses from the work of [START_REF] Moretti | High Strain Rate Deformation Behavior and Recrystallization of Alloy 718[END_REF] for the annealed specimens using the standard crystal plasticity model are shown in Fig. 16a for three temperatures. The material constants and fitted material parameters are summarized in Table 4, after [START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF]. In addition, the material parameters τ 0 , K, m from ( 12) are considered as temperature-dependent. These calibrated temperature-dependent material parameters are given in Table 5. Moreover, the constant parameter d c in the dislocation density-based model, which controls the saturation of the stress-strain behavior, is also taken as temperature-dependent (Table 4). The initial dislocation density ρ r 0 (= ϱ r 0 /b 2 ) is assumed to be 10 11 m -2 and chosen the same for all slip systems. As a general trend, with increasing temperature, the flow strength of Inconel 718 and consequently the plastic work and dislocation density decrease. In Fig. 16b the dislocation density rapidly increases at an initial deformation stage and saturates during further straining. Moreover, samples deformed at lower temperatures show higher dislocation densities than those deformed at higher temperatures. The material parameters of Inconel 718 are calibrated for the aged samples against the experimental work of [START_REF] Iturbe | Mechanical characterization and modelling of inconel 718 material behavior for machining process assessment[END_REF]. The stored energy evolution displayed in Fig. 17a shows a similar trend as the dislocation density. It increases rapidly at the initial stage of deformation and saturates with a further increase in deformation. The evolution of TQCs is provided in Fig. 17b. The predicted lowest value of TQC using ξ = 10 is 0.825. In comparison, in the previous section, polycrystalline austenitic steel shows a lowest TQC value of 0.58 with ξ = 10 (Fig. 15), as calibrated from experimental measurements. Unfortunately 17) κc (fitted) Eq. ( 16) dc (fitted) (293K) Eq. ( 16) dc (fitted) (673K) Eq. ( 16) dc (fitted) (1073K) Eq. ( 16) h0 Eq. ( 15 Fig. 16: Annealed Inconel 718 under compressive loading: (a) Predicted averaged stress-strain curves at high strain rates and three initial temperatures (293 K-1073 K) and comparison against the experimental data from [START_REF] Moretti | High Strain Rate Deformation Behavior and Recrystallization of Alloy 718[END_REF] and (b) predicted evolution of the averaged dislocation density (Eq. ( 58)) as a function of strain.

no such calorimetric measurements are available for Inconel 718 so that no definitive conclusion can be drawn from the predicted values of TQC.

Contribution of GNDs to stored energy

In this section, the stored energy and TQC are predicted using the micromorphic crystal plasticity model described in section 2.5. The model includes additional stored energy due to the norm of microslip gradient, see [START_REF] Maugin | The Thermomechanics of Nonlinear Irreversible Behaviors[END_REF]. The contribution K • A ∼ • K /2 = A∥K ∥ 2 /2 involving the higher order modulus A, is an estimation of GND density through the norm of the slip gradient vector grad γ χ ≃ grad γ cum , according to [START_REF] Ling | A reduced micromorphic single crystal plasticity model at finite deformations. application to strain localization and void growth in ductile metals[END_REF]. The micromorphic model not only includes an evaluation of the GND density but also a size-dependent additional hardening term in [START_REF] Shen | Stored energy density research on the fatigue crack initiation at twin boundary and life prediction of Inconel718 superalloy[END_REF]. This additional 56), ( 57)) and (b) Taylor-Quinney coefficient (Eq. ( 61), ( 62)) as functions of strain for Inconel 718 under compressive loading using ξ = 10.

hardening term -Div M = -A∆γ χ is proportional to the Laplacian of microslip distribution, in a similar way as in the Aifantis strain gradient plasticity model [START_REF] Aifantis | The physics of plastic deformation[END_REF].

The analytical expressions used for the prediction of stored energy density and TQC are given in Eq. ( 63) and [START_REF] Gérard | Modeling of latent hardening produced by complex loading paths in fcc alloys[END_REF], respectively. They are averaged over the entire polycrystalline volume. The stored energy evolution predicted by the micromorphic model is compared with the results obtained previously with the standard crystal plasticity model.

The material constants and calibrated material parameters of the standard crystal plasticity part of the model are taken from Table 3 for austenitic steel 316L. The micromorphic part of the model involves two additional parameters: A and H χ . A high value H χ = 50000 MPa of the penalty parameter is taken constraining the microslip variable γ χ to almost coincide with the cumulative slip variable γ cum . Under these conditions, the micromorphic model can be regarded as a strict strain gradient plasticity model as discussed in [START_REF] Ling | A reduced micromorphic single crystal plasticity model at finite deformations. application to strain localization and void growth in ductile metals[END_REF]. A precise identification of the micromorphic model requires first the calibration of the standard crystal plasticity model from single crystal material behavior. Then, the gradient parameter A is derived from grain size dependent tensile curves [START_REF] Cordero | Generalised continuum modelling of grain size effects in polycrystals[END_REF]. This strategy has not been undertaken for the considered material due to limited available data for single crystal austenitic steel 316L. Instead, the gradient parameter A is chosen as 0.02 N such that predicted stress-strain response remains close to experimental response obtained from the work of [START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF] keeping the standard crystal plasticity model parameters previously identified for a polycrystal with a given grain size of 80 µm.

The fitted average stress-strain response against the experimental data is shown in Fig. 18a. A characteristic length scale of the model can be defined from the material parameters as ℓ =

A(H+Hχ) |H|Hχ

where H is the plastic hardening modulus. This is a varying length scale as the hardening modulus H is not constant [START_REF] Ling | A reduced micromorphic single crystal plasticity model at finite deformations. application to strain localization and void growth in ductile metals[END_REF]. The calculated characteristic length scale is about 2 µm for A = 2 × 10 -2 N, H χ = 5 × 10 4 MPa and H = 3000 MPa.

Grain boundary conditions play an essential role in the micromorphic polycrystal response. These interface conditions are chosen as intermediate between so-called microfree and microhard conditions, as done in [START_REF] Forest | Cosserat modelling of size effects in the mechanical behaviour of polycrystals and multiphase materials[END_REF][START_REF] Cordero | Grain size effects on plastic strain and dislocation density tensor fields in metal polycrystals[END_REF]. The microslip variable γ χ is assumed to be continuous at the interface while surface traction T and generalized surface traction M in Eq. ( 20) are also continuous. In contrast, microfree or microhard conditions respectively assume vanishing generalized traction vector or microslip along the grain boundaries. They are known to predict lower and upper bounds of the size-dependent material response, respectively [START_REF] Bargmann | Modeling of polycrystals with gradient crystal plasticity: A comparison of strategies[END_REF]. Intermediate material responses are obtained using the present continuity conditions.

Stored energy density is predicted using an ad-hoc factor ξ of 1. The comparison of the predicted stored energy using the standard and micromorphic crystal plasticity 63), ( 64)), and SSDs only (Eq. ( 56), ( 57)).

The stress-strain responses are validated against the experimental results from [START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF] for austenitic steel 316L. models is shown in Fig. 18b. It shows that the predicted volume averaged stored energy using the micromorphic crystal plasticity model is higher than that of the standard crystal plasticity model. The standard crystal plasticity model lacks a characteristic length scale that is associated with the GNDs. On the other hand, the presence of characteristic length scale in the micromorphic crystal plasticity model takes into account the contribution of GNDs along with SSDs and therefore gives an increased total dislocation density and, consequently, stored energy. With the chosen material parameters, the GND contribution is no more than one third of the SSD one. The stored energy expression ( 63) is composed of three contributions in the form:

E s = 1 2 µξ N r=1 (ζ r ) 2 Ec + 1 2 H χ e 2 p Ep + 1 2 AK • K Eg . ( 67
)
These contributions are illustrated in Fig. 19a. It is observed that E c avg based on scalar dislocation densities, is larger than E g avg . The contribution E p avg is specific to the micromorphic model and is found to be negligible (in this respect, see also [START_REF] Abatour | Toward robust scalar-based gradient plasticity modeling and simulation at finite deformations[END_REF] in the isotropic plasticity case). Finally, E c avg > E g avg > E p avg . Additional simulations are performed with microhard grain boundary conditions. The microhard grain boundary condition corresponds to vanishing microslip γ χ at the grain boundaries, i.e. γ χ = 0. These grain boundary conditions are prescribed by setting the nodal values of γ χ on the grain boundary surfaces to zero. The comparison of the obtained responses using microhard conditions with higher order moduli A = 2 × 10 -2 N and 2 × 10 -5 N to that of the results from the experimental test is presented in Fig. 19b. The response obtained using A = 0.02 N is too strong compared to experimental results. On the other hand, A = 2 × 10 -5 N gives a significantly lower response than A = 2 × 10 -2 N. This exaggerated effect of grain boundary conditions is probably due to the fact that the material parameters of the standard part of the micromorphic crystal plasticity model were not identified from single crystal data for austenitic steel but from polycrystal data at a given grain size. This means that the parameters of the standard crystal plasticity model already contain some component of grain size hardening. A comparison of the predicted total stored energy with intermediate and microhard grain boundary conditions is drawn in Fig. 20. The Stored energy (MJm -3 )

Microhard(A = 2x10 -2 N ) Microhard(A = 2x10 -5 N ) Intermediate(A = 2x10 -2 N )
Fig. 20: Comparison of the predicted volume averaged stored energy (Eq. ( 63), ( 64)) using microhard and intermediate grain boundary conditions with ξ = 1.

predicted total stored energy with the microhard grain boundary conditions is significantly higher than that of the intermediate grain boundary conditions. These values lead however to unrealistic stress levels, as seen in Fig. 19b. As a conclusion, the crystal plasticity model identified from the polycrystal experimental response, predicts low storage of energy, even after adding the GND content predicted by the strain gradient plasticity approach. There remain possibilities of increasing the amount of stored energy by proper identification of the crystal plasticity model from single crystal data.

Conclusions

The numerical simulation work presented is intended to provide a first insight into the prediction of the stored energy and evolution of TQC in single and poly-crystalline FCC metallic materials based on thermodynamically consistent formulations of standard and micromorphic crystal plasticity models. The following conclusions are drawn from the study:

• Thermodynamic consistency requires checking the positivity of the entropy production rate predicted by the model. For that purpose, it is necessary to select internal variables and define the Helmholtz free energy potential, which is rarely done in the literature. A phenomenological model was described having the advantage that positive dissipation rate for all loading conditions can be enforced by introducing a suitable dissipation potential. In contrast, the introduction of physically based internal variables with more complex evolution equations is not compatible with the existence of a dissipation potential in general. We analytically demonstrated using simplified cases that the positivity of the dissipation rate in dislocation densitybased models can be ensured by imposing constraints on the model parameters.

In other cases, positivity of the dissipation rate must be checked at each material point and at each time increment, as done in the presented simulations. • The stored energy and TQC for aluminum and copper single crystals (<001> and <111> crystal orientations) has been predicted using a standard crystal plasticity model. The found values are higher than usual values close to 0.8 -0.9 and than the experimentally measured average TQC values from the work of [START_REF] Rittel | On the dynamically stored energy of cold work in pure single crystal and polycrystalline copper[END_REF] for copper single crystals. • Polycrystal simulations with a value ξ = 1 underestimates the predicted stored energy compared to the experimentally measured values for austenitic steel. On the other hand, predicted stored energy with ξ = 10 shows good agreement with the experimentally measured values, which is consistent with the single crystal identification. • The contribution of GNDs was put forward as a part of the studied dark energy, within a thermodynamically consistent formulation of micromorphic/SGP crystal plasticity. Based on the presented identification strategy of the SGP material parameters, it is found that the predicted stored energy is only moderately higher than the standard estimation. • Grain boundary conditions play an important role in SGP simulations. The microhard grain boundary conditions induce too strong stress response compared to experimental results. The higher order modulus A was calibrated based on intermediate interface conditions enforcing continuity of slip. There is hope that a higher content of stored energy could be gained by a more accurate identification of the micromorphic model from single crystal data combined with grain size dependent experimental results. However enriching the GND content will be at the expense of the SSD part and the increase in stored energy is expected to remain limited, still leaving dark energy contributions to be explored.

As such, this work leaves more questions than it answers. A primary objective is to show that frequently used dislocation based crystal plasticity models can be made thermodynamically consistent so that stored energy can indeed be predicted. The results indicate that the storage of energy predicted by this class of models is underestimated. The high values of the ad-hoc parameter ξ show the weakness of the proposed thermodynamic crystal plasticity framework and ask for explanations for the origin of the dark energy. This may indicate that the dislocation density, an average quantity for SSD and GND, is not sufficient to evaluate stored energy. Additional sources of concentrated energy in the form of specific dislocation structures exist and should be incorporated in the model. This calls for more accurate estimates of energy stored by dislocations. Some improvements already exist in the literature. For instance, the authors in [START_REF] Saimoto | Forensic analyses of microstructure evolution of stage II & III: new assimilated model for work-hardening in FCC metals[END_REF] insist on the importance of a detailed description of coplanar annihilation of dislocations. Regarding GND contributions, a more precise form of the energy in stacked pile-ups is studied in [START_REF] Mesarovic | Thermodynamic coarsening of dislocation mechanics and the size-dependent continuum crystal plasticity[END_REF][START_REF] Baskaran | Energies and distributions of dislocations in stacked pile-ups[END_REF]. Alternative more elaborate thermodynamic settings have been proposed such as the two temperature model in [START_REF] Lieou | Thermomechanical conversion in metals: dislocation plasticity model evaluation of the Taylor-Quinney coefficient[END_REF]. These should be explored and compared to the present standard formulations. Where is the dark energy hidden in dislocation microstructures? Local interactions between dislocations as computed by discrete dislocation dynamics (DDD) approaches [START_REF] Zbib | A multiscale model of plasticity based on discrete dislocation dynamics[END_REF] may well play a significant role and cannot be captured by functions of the sole total dislocation densities. Luo et. al. [START_REF] Li | Stable stacking faults bounded by Frank partial dislocations in Al7075 formed through precipitate and dislocation interactions[END_REF] recently proposed that vacancy formation as a result of interaction between moving dislocations (induced by deformation) and precipitates can be a source of energy storage high enough for stacking fault formation. Stacking fault energy can be regarded as an additional source of stored energy storage. It is also worth noting that, according to the results from Rittel et al. [START_REF] Rittel | The dependence of the Taylor-Quinney coefficient on the dynamic loading mode[END_REF], the TQCs of pure titanium (except for tension test) and 1020 steel are about 0.8 -0.9, similar to the theoretical value of 0.9, unlike the other precipitation hardening alloys, Al2024 and Al5086, with TQCs around 0.2.

It should be noted also that the micromorphic model used involves only the gradient of the cumulative slip variable. Consideration of the full dislocation density tensor may help improving the predictions [START_REF] Phalke | Modeling size effects in microwire torsion: A comparison between a lagrange multiplier-based and a CurlF p gradient crystal plasticity model[END_REF]. Grain boundaries may also be places where dark energy is hidden. More elaborate treatment of grain boundary conditions than the presented intermediate or microhard contitions could help capture these contributions. Finite resistance against dislocation gliding is introduced in [START_REF] Fredriksson | Competition between interface and bulk dominated plastic deformation in strain gradient plasticity[END_REF][START_REF] Ekh | Influence of grain boundary conditions on modeling of size-dependence in polycrystals[END_REF][START_REF] Husser | Size affected dislocation activity in crystals: Advanced surface and grain boundary conditions[END_REF] and gives more realistic dislocation-grain boundary interactions. More general grain boundary conditions have been proposed in [START_REF] Wulfinghoff | Equivalent plastic strain gradient crystal plasticityenhanced power law subroutine[END_REF] allowing for the transition from microhard conditions to microfree (or constant generalized tractions) once a threshold is reached at the grain boundary. The contribution of grain boundary energy itself, as introduced in [START_REF] Aifantis | The role of interfaces in enhancing the yield strength of composites and polycrystals[END_REF], is significant for nano-grains and becomes negligible for conventional engineering grains sizes. The strain gradient plasticity model used in this work accounts for size effects and can therefore be applied to study the effect of grain size on energy storage in austenitic steel, as done experimentally in [START_REF] Oliferuk | Effect of the grain size on the rate of energy storage during the tensile deformation of an austenitic steel[END_REF]. However, the current identification of the strain gradient plasticity model is too crude to make reliable predictions. We leave for future work a more systematic analysis of grain size effects. The question of the distribution of stored energy inside the grains is even more challenging although first intragranular thermal measurements are available such as [START_REF] Charkaluk | Microplasticity in Polycrystals: A Thermomechanical Experimental Perspective[END_REF]. Finally, anisotropic crystallographic textures, strain localization phenomena and recrystallization were not considered in the analysis and are known to play a significant role in energy storage.

U 2 = 0 (U 2 ( 3 e 2 e 1 Fig. 1 :

 2023211 Bottom surface) Top surface)U 1 = 0 (Back surface) U 3 = 0 (Side surface)e Applied boundary conditions on cube volume element for single and polycrystalline simulations.

Fig. 2 :Fig. 3 :

 23 Fig. 2:Comparison of stress-strain responses for <001> and <111> crystal orientations at a temperature of 298 K: (a) Aluminum single crystals validated against the experimental results of[START_REF] Hosford | Tensile deformation of aluminum single crystals at low temperatures[END_REF], and (b) copper single crystals validated against the experimental results of[START_REF] Takeuchi | Work hardening of copper single crystals with multiple glide orientations[END_REF].

Fig. 6 :

 6 Fig. 6: (a) Comparison of stress-strain responses against the experimental data obtained by [8] for single crystal copper (<123> orientation) under compression loading using boundary conditions of Fig. 1; (b) Predicted dislocation density ρ r (= ϱ r /b 2 ) evolution at two different strain rates.

Fig. 8 :Fig. 9 :

 89 Fig. 8: Computing polycrystals with structured FE meshes containing 64 cuboidal copper grains with random orientations: (a) 1 element per grain, (b) 8 elements per grain and (c) 27 elements per grain. Colors represent the individual grains.

Fig. 10 :

 10 Fig. 10: Computing polycrystals with FE meshes of 64 randomly oriented grains generated by Voronoï tessellation: (a) Coarse (1774 nodes, 8274 elements) and (b) fine (16402 nodes, 85251 elements) meshed. Colors represent the individual grains.

Fig. 11 :Fig. 12 :

 1112 Fig. 11: RVE generated by Voronoï tessellation: (a) 64 grains, (b) 125 grains, and (c) 200 grains. Color represents the individual grains.

Fig. 13 :

 13 Fig.13: Volume averaged stress-strain responses obtained using copper grains created by Voronoï tessellation to study the effect of (a) mesh size (RVE of 64 grains) and (b) total number of grains in the RVE.

Fig. 14 :

 14 Fig.14: Polycrystalline austenitic steel 316L under tensile loading at room temperature: (a) Averaged stress-strain curve validated against the experimental work of[START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF], and (b) predicted volume averaged stored energy over whole FE model (ξ = 1 and 10) (Eq. (56), (57)) and comparison against the experimental measurements from the work of[START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF].

10 Fig. 15 :

 1015 Fig. 15: Predicted the volume averaged Taylor-Quinney coefficient (Eq. (61), (62)) using numerical simulations with two different values of an ad-hoc factor (1 and 10) for austenitic steel 316L.

, 1900 s - 1 Simulation, 293K, 1900 s - 1 Experiment, 673K, 1900 s - 1 Simulation, 673K, 1900 s - 1 Experiment, 1073K, 2000 s - 1 Simulation 14 293K, 1900 s - 1 673K, 1900 s - 1 1073K, 2000 s - 1 (

 1111114111

293K, 1900 s - 1 673K, 1900 s - 1 1073K, 2000 s - 1 (Fig. 17 :

 11117 Fig. 17: Predicted evolution of volume averaged (a) stored energy (Eq. (56), (57)) and (b) Taylor-Quinney coefficient (Eq. (61), (62)) as functions of strain for Inconel 718 under compressive loading using ξ = 10.

3 )Fig. 18 :

 318 Fig. 18: Polycrystalline austenitic steel under tensile loading at room temperature: (a) Averaged stress-strain responses predicted considering the contribution of SSDs and GNDs (micromorphic model), and SSDs only (standard crystal plasticity model); (b) Predicted evolution of volume averaged stored energy (using ξ = 1) considering the contribution of SSDs and GNDs (Eq. (63), (64)), and SSDs only (Eq. (56), (57)). The stress-strain responses are validated against the experimental results from[START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF] for austenitic steel 316L.

Fig. 19 :

 19 Fig. 19: (a) Contribution of each term to the total stored energy in micromorphic crystal plasticity model (Eq. (67)); (b) Predicted stress-strain responses using microhard grain boundary condition using A = 2×10 -2 N and 2×10 -5 N. The experimental stress-strain responses are from taken from[START_REF] Oliferuk | Rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel[END_REF] for austenitic steel 316L.

Table 2 :

 2 Material constants and fitted material parameters for the simulation of single crystalline copper under tensile loading. The interaction matrix coefficients (h 0 -h 5 ) are taken from[START_REF] Kubin | Modeling dislocation storage rates and mean free paths in face-centered cubic crystals[END_REF].

			1/m	5	27
	h0	h1	h2	h3	h4	h5	ξ
	Eq. (15)	Eq. (15)	Eq. (15)	Eq. (15)	Eq. (15)	Eq. (15)	Eq. (56)
	0.122	0.122	0.07	0.625	0.137	0.122	1
	C11	C12	C44	ρ Eq. (60)

Table 4 :

 4 Numerical values of material parameters used for the numerical simulation of polycrystalline Inconel 718 under compressive loading, after[START_REF] Phalke | Adiabatic shear banding in fcc metallic single and poly-crystals using a micromorphic crystal plasticity approach[END_REF].

	C11	C12	C44	ρ Eq. (60)	Cε Eq. (60)	b	µ
	194 GPa	142 GPa	90 GPa	7800 kgm -3	435 Jkg -1 K -1	0.249 nm 77.2 GPa
	ϱ r 0 (fitted)
	Eq. (

Table 5 :

 5 Temperature dependent material parameters used in the numerical simulations for annealed Inconel 718.

			Temperature (K)	τ0 (MPa) Eq. (15)	K(MPa.s 1/m ) Eq. (12)	m Eq. (12)
			298		210	5	10
			923		150	5	10
			1073		80	8	6
			1173		60
			1323		10	40	4
			1523		1	50	2
	0.00	0.05	0.10	0.15	0.20

Acknowledgement

The research leading to these results has received funding from the European Union's Horizon 2020 research and innovation programme, Marie Sklodowska-Curie Actions,under grant agreement no. 764979, project "ENABLE", "European Network for Alloys Behaviour Law Enhancement".

Data and code availability

All data are contained in the paper or accessible upon request. The finite element code Zset is available by transvalor.com

Declarations

Author contributions

VP contributed to conceptualization, investigation, data curation, formal analysis, and writing-original draft; SF contributed to investigation, methodology, and writingoriginal draft; HJC contributed to formal analysis, and writing-original draft; TR contributed to investigation, resources, and writing-review and editing; AR contributed to investigation, resources and writing-review and editing. All authors have read and agreed to the published version of the manuscript.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper

Supplementary information

Not applicable.

Ethical approval

Not applicable.