Y F Zhang

R Vicquelin
email: ronan.vicquelin@ecp.fr

Controlling bulk Reynolds number and bulk temperature in channel flow simulations

Keywords: Channel flow, Control, Bulk Reynolds number, Bulk temperature

Bulk Reynolds number and bulk temperature are key quantities when reporting results in channel flow simulations. There are situations are when one wishes to accurately control these parameters while changing some numerical or physical conditions. A method to control the bulk Reynolds number and the bulk temperature in channel flow simulations is detailed. An ordinary differential equation is prescribed for the additional source term in the momentum balance equation so that the transient regime of the simulation is thoroughly tuned in order to efficiently and accurately reach the target Reynolds number value. A similar treatment is applied for the additional volume heat source term in the energy balance equation. The proposed method is specifically interesting when studying complex multi-physics in channel flow configurations when non-dimensionalization of the equations is no longer practical.

Introduction

Channel flow is a simple configuration that has been widely studied numerically in order to analyze turbulent boundary layers and validate models [START_REF] Moin | Direct numerical simulation: A tool in turbulence research[END_REF][START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF]. For fully developed channel flows, periodic conditions are considered in the infinite spanwise direction, and also along the streamwise direction. A homogeneous source term S i is then added to the momentum equations to compensate for viscous forces and drive the flow at a given bulk velocity. In direct numerical simulations, the momentum equation in the streamwise direction (i = 1) reads

∂ ∂t (ρu 1 ) + ∂ ∂x j (ρu j u 1 ) = - ∂p ∂x 1 + ∂τ 1 j ∂x j + S 1 (1)
where S 1 is the homogeneous source term in the streamwise direction while S 2 = S 3 = 0. The source term S 1 enforces the flow mass flow rate and then determines the obtained Reynolds number. In standard flows with constant flow properties and without multiple physical phenomena (chemistry, radiation, ...), non-dimensional equations can be written and this source term is directly related to the intended bulk or friction Reynolds number. Alternatively, the source term can be simply be updated step by step to compensate for the change of the intended mass flux [START_REF] Garcia-Mayoral | Hydrodynamic stability and breakdown of the viscous regime over riblets[END_REF]. However, in more complex flows featuring variable properties due to an explicit dependency on temperature for example [START_REF] Bocquet | A compressible wall model for large-eddy simulation with application to prediction of aerothermal quantities[END_REF], or involving multiphysics such as chemistry [START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF] or radiation [START_REF] Ghosh | Effects of radiative heat transfer on the structure of turbulent supersonic channel flow[END_REF][START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF], equations are kept in their dimensional form and S 1 must be determined differently to reach a target bulk Reynolds number Re t b . A first method to determine S 1 consists in taking a fixed constant in time S ref [START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF][START_REF] Ghosh | Effects of radiative heat transfer on the structure of turbulent supersonic channel flow[END_REF] that is either chosen arbitrarily or more carefully evaluated from friction coefficient formulae that depend on the Reynolds number. Usually, the considered formulae are accurate for simple flows (constant properties, no multiphysics) but can become considerably erroneous as the studied flow is more and more complex. In such a case, the final Reynolds number that is obtained is different from the intended one.

A second kind of method consists in dynamically adapting the source term value after each iteration so that the Reynolds number is brought towards its target value [START_REF] Lenormand | Large eddy simulation of subsonic and supersonic channel flow at moderate Reynolds number[END_REF][START_REF] Bocquet | A compressible wall model for large-eddy simulation with application to prediction of aerothermal quantities[END_REF][START_REF] Cabrit | Direct simulations for wall modeling of multicomponent reacting compressible turbulent flows[END_REF][START_REF] Malaspinas | Wall model for large-eddy simulation based on the lattice boltzmann method[END_REF]. However, none of these methods have been carefully characterized with a dedicated study. The procedures reported in [START_REF] Bocquet | A compressible wall model for large-eddy simulation with application to prediction of aerothermal quantities[END_REF][START_REF] Malaspinas | Wall model for large-eddy simulation based on the lattice boltzmann method[END_REF] are similar to

S n+1 1 = S ref + ρ t b u t b -V ρ n u n dV V τ ref , (2)
where S n+1 1 is the updated value of the uniform momentum source term at iteration n + 1, ρ t b is the target value of bulk density and u t b is the target value of bulk velocity. Equation 2 describes the adaptation of S n+1 1 because of the difference between the bulk mass flux at iteration n and its target value ρ t b u t b . This difference is related to the difference between the simulation bulk Reynolds number and the target value Re t b . The relaxation time τ ref used in Eq. 2 is typically expressed in terms of the channel time scale δ/u τ . When using Eq. 2, the permanent regime is reached after a transient stage and S n+1 1 tends towards a constant value denoted by S 1 , and

ρ t b u t b -V ρ n u n dV V = S 1 -S ref τ ref . (3)
Since the obtained stationary value S 1 is different from the empirically determined S ref in complex flows, a finite bias between Re b and Re t b unfortunately remains in such flows. This biais is not present in the formulation from Lenormand et al. [START_REF] Lenormand | Large eddy simulation of subsonic and supersonic channel flow at moderate Reynolds number[END_REF] which can be reformulated as a PI controller with constant coefficients and a time response close to the computation time step. The efficiency of the approach was established a posteriori.

Equation 2 was thus a first attempt to control the bulk Reynolds number. This objective is here pursued by proposing an approach which ensures that the channel flow converges exactly and efficiently to the target bulk Reynolds number a priori and a posteriori. The approach appears as a modified PI controller with time varying coefficients. Carefully tuning the time varying coefficients allows to derive a second order ordinary differential equation with constant coefficients for Re b whose time response can then be exactly controlled. The method is described in the following section. Its efficiency is then demonstrated in several channel flow configurations. Similarly, the method is derived and applied for the equivalent control of bulk temperature when studying turbulent heat transfer. Finally, a strategy for the control of both bulk Reynolds number and bulk temperature with variable thermo-physical properties accounted for is derived and validated in the section 4.

Control of the bulk Reynolds number

Formulation

Integration of Eq. 1 over the whole computational domain V gives

d dt V ρu 1 dV = τ 1 j n j dS + S 1 V (4)
where n j is the outward surface normal vector. The integration of the pressure gradient and the convective terms is null because of the applied periodic boundary conditions in the X and Z directions.

The integrated term on the left side of Eq. 4 is related to the bulk Reynolds number Re

b = ρ b u b δ µ b , V ρu 1 dV = ρ b u b V = µ b Re b V δ ( 5
)
where δ is the half-width of the channel and the subscript b is related to the aforementioned bulk quantities. The integrated wall shear stress is split into two contributions,

τ 1 j n j dS = -(τ w,1 + τ w,2 )S w (6)
as τ w,1 and τ w,2 , the average wall shear stresses on the lower and upper wall respectively, can be different in the general case. S w denotes the surface area of each wall.

Neglecting variations of the bulk dynamic viscosity and noticing that V = 2δS w , Eq. 4 becomes

dRe b dt = - τ w,1 + τ w,2 2 µ b + δ µ b S 1 . (7)
Differentiating Eq. 7 to make dS 1 dt appear gives

d 2 Re b dt 2 = - d dt τ w,1 + τ w,2 2 µ b + δ µ b dS 1 dt . (8)
As in many control algorithms, we choose the evolution of the homogeneous source term S 1 to be determined by the dynamics of the controlled quantity Re b :

dS 1 dt = a 1 (t) dRe b dt + a 2 (t), (9)
where a 1 (t) and a 2 (t) are to be determined. Equation 8is then written as

d 2 Re b dt 2 -a 1 (t) δ µ b dRe b dt -a 2 (t) δ µ b + d dt |τ w,1 | + |τ w,2 | 2 µ b = 0. ( 10
)
In order to obtain a controlled relaxation towards the target value Re t b , Eq. 10 is transformed into the following second-order ordinary differential equation with constant coefficients α and β for the difference (Re b -Re t b ),

d 2 dt 2 (Re b -Re t b ) + α d dt (Re b -Re t b ) + β (Re b -Re t b ) = 0, (11)
by setting a 1 (t) and a 2 (t) as

a 1 (t) = -α µ b δ , (12)
a 2 (t) = µ b δ β (Re t b -Re b (t)) + d dt |τ w,1 | + |τ w,2 | 2 µ b . ( 13
)
The dynamics of Re b (t) is then controlled by the discriminant ∆ = α 2 -4β. Among the different possible regimes, the critically damped regime (∆ = 0) gives the shortest decay time in order for the controlled bulk Reynolds number to reach the desired targeted value. It is therefore chosen, leading to α = 2β 1/2 . Defining τ = 2 α introduces the controlling decay time that can be chosen arbitrarily. Accounting for the initial conditions Re b t=0 and dRe b dt t=0 , the temporal evolution obtained for the bulk Reynolds number is

Re b (t) = Re t b + Re b t=0 -Re t b 1 + t τ + dRe b dt t=0 t e -t τ , (14)
The final formulation for the source term is therefore:

dS 1 dt = - µ b δ 2 τ dRe b dt + Re b (t) -Re t b τ 2 - d dt τ w,1 + τ w,2 2 µ b , (15)
which appears as a modified PI controller with thoroughly tuned coefficients. Indeed, equation 15 contains time varying terms (apart from Re b (t) and dRe b dt ) instead of constants. Beside, their values are determined from the exact evolution law of the bulk Reynolds number rather than by linearization of the system equations around a mean state. Hence, any non-linear effects in the physical response due to a change in the source term is accounted for. The temporal evolution described by Eq. 14 is exact when bulk dynamic viscosity is constant. This is the only assumption that has been made. Even when this is not true, Eq. 15 ensures that the statistically steady state will always correspond to Re b = Re t b . Initial condition S 1 (0) could be determined from Eq. 7 to impose an smooth departure from the initial Reynolds number by imposing dRe b dt t=0 = 0. Alternatively, the initial source term value is here fixed to

S ref = τ ref w
δ a reference value that estimates the source term for the targeted Reynolds number. Equation 15 is coupled with the set of Navier-Stokes conservative equations. Using forward Euler time integration, the value of the momentum source term is updated at each new iteration. This approach has been applied successfully in direct numerical simulations and large-eddy simulations of channel flows with or without radiative energy transfer for identical bulk Reynolds number [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF][START_REF] Zhang | A wall model for les accounting for radiation effects[END_REF].

Sometimes, the source term in the momentum equation is expressed as a body force ρS 1 . In variable density flows, the source term ρS 1 is then not uniform but the proposed approach to determine S 1 can be simply adapted to this case.

Validation

The computational domain is shown in figure 1 is 2πδ × 2δ × πδ. The mesh resolution is the same as in case C3 reported in [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF]. The set of governing equations is solved with the finite-volume solver YALES2 [START_REF] Moureau | Design of a massively parallel CFD code for complex geometries[END_REF][START_REF] Moureau | From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: Filtered laminar flame-pdf modeling[END_REF] under a low Mach-number approximation. As detailed in [START_REF] Zhang | Physical study of radiation effects on the boundary layer structure in a turbulent channel flow[END_REF], the numerical setup is composed of a centered fourth-order spatial discretization and a fourth-order time integration.

Isothermal case

The only assumption made in the previous derivation is neglecting the variations of the bulk dynamic viscosity. The first studied case is isothermal (T = 1000 K) and the constant viscosity assumption is correct. Three different decay times τ are considered: 0.33δ/u re f τ , 0.033δ/u re f τ , and 0.0033δ/u re f τ . They are set as proportionate to the estimated channel time scale δ/u re f τ . The initial state of the channel flow simulation is a statistically steady turbulent flow at Re b = 5850. The control method is applied to obtain Re b = Re t b = 7000. The mesh resolution remains adequate for this larger Reynolds number.

Figure 2 (a) shows that, whatever the imposed time scale τ, the temporal evolution of the bulk Reynolds number expressed in non-dimensional time units t/τ follows the theoretical solution given by equation 14. It is then possible to reach the target value Re t b at any speed. Obviously, the time scale must still be chosen in order not to penalize the stability of time integration of the full set of equations. The very small deviation from the theoretical solution of the case with the smallest time response is due to numerical errors in the time integration because of the time decay τ becoming of the same order as the simulation time step.

Figure 2 (b) shows the corresponding temporal evolution of the friction Reynolds number defined as Re τ = ρ w u τ δ/µ w where u τ is the actual friction velocity averaged in homogeneous directions on both walls, ρ w and µ w are the density and dynamic viscosity at the wall. Re τ is here plotted against the normalized time t u re f τ /δ to highlight that, for short decay times τ, the adaptation of the flow field to the change in the bulk Reynolds number takes a given amount of time that roughly scales with u τ /δ.

Case with variable flow properties

Another computational case is now considered to study the impact of the constant bulk viscosity assumption made to derive the controlling method. Wall temperatures are set to T w,c = 950 K and T w,h = 2050 K for the cold and hot wall, respectively. From the instantaneous turbulent field of temperature, the dynamic viscosity µ is computed from a tabulated function of temperature built from detailed computation with the CHEMKIN package [START_REF] Kee | A fortran computer code package for the evaluation of gas-phase, multicomponent transport properties[END_REF][START_REF] Kee | CHEMKIN-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics[END_REF]. This time, the initial state for the control method is a turbulent channel flow at Re b = 5850 whose statistically steady state is reached with the modified wall temperatures T w,c and T w,h . The target value for the bulk Reynolds number is set to Re t b = 7000. The equation for dS 1 dt can be adapted to account rigorously for a variable bulk viscosity (see section 4). It was chosen here to rely instead on the robustness of the derived controller. Results are shown in figure 3 for two different decay times τ. In this case, the evolution of the bulk Reynolds number follows the ideal response with very small deviations. Then, the target value of the bulk Reynolds number is effectively reached as imposed by Eq. 15. Because of the different density and viscosity at each wall, a friction Reynolds number Re τ is defined for each side. As observed in the isothermal case, the adaptation of the mean wall shear stress takes a given physical time different from the prescribed value of τ.

Bulk Temperature

The analogical issue of bulk Reynolds number control raises when studying compressible flows or turbulent heat transfer in channel flow simulations. The bulk temperature or enthalpy is then controlled by the wall heat transfer. When both walls have a different temperature or when heat is generated by viscous dissipation in high-speed flows, the bulk temperature is physically determined by the problem setup. In other configurations ( [START_REF] Coleman | A numerical study of turbulent supersonic isothermal-wall channel flow[END_REF] for example), one might want to control the bulk temperature by adding a new source term in the enthalpy or energy balance equation. The introduction of an energy source term is also required to maintain the bulk temperature T b in periodic simulations when considering identical temperatures at both walls that are different from T b without any physical heat source term. By adding an energy source term Ω, the enthalpy transport equation is given by

∂ ∂t (ρh) + ∂ ∂x j (ρu j h) = - ∂q j ∂x j + ∂p ∂t + u j ∂p ∂x j + τ i j ∂u i ∂x j + Q + S i u i Qtot +Ω ( 16
)
where pressure effects, viscous heating, additional physical energy source term Q and the work S i u i due to the added momentum source term are considered in the general case and gathered in Qtot . The previous method is now derived for the controlling source term Ω. It is first validated when temperature is considered as a passive scalar field and where the bulk Reynolds number is fixed. Then, results for simultaneous control of the bulk Reynolds number and temperature are presented.

Formulation

The control method is obtained similarly to the one developed for the bulk Reynolds number, i.e. it is exact for a flow with constant properties. More complex expressions are given in section 4 by accounting for effects of variables flow properties. Integrating the enthalpy balance equation gives

d dt V ρhdV = -q j n j dS + V Qtot dV + ΩV, ( 17
)
where a bulk enthalpy h b can be introduced, with

V ρhdV = ρ b h b V. ( 18
)
When variable properties are considered, the thermal capacity at constant pressure c p is not constant and the bulk temperature T b can be defined by inverting the relation h b = h(T b ). Neglecting temporal variations of bulk density, equation 17 becomes

ρ b V dh b dt = -q j n j dS + V Qtot dV + ΩV. (19)
The integrated wall heat flux is split into two parts,

q j n j dS = (q w,1 + q w,2 )S w (20)
where q w,1 and q w,2 are the average wall conductive fluxes on the lower and upper wall, respectively. Neglecting variations of c p gives, dh b dt = c p b dT b dt , and the temporal evolution of the bulk temperature is described by dT

b dt = - q w,1 + q w,2 2δ ρ b c p b + 1 ρ b c p b V V Qtot dV + 1 ρ b c p b Ω (21)
Differentiating in respect to time gives Taking the following form for the source term Ω,

d 2 T b dt 2 = - 1 2δ d dt q w,1 + q w,2 ρ b c p b + d dt 1 ρ b c p b V V Qtot dV + 1 ρ b c p b dΩ dt (22)
dΩ dt = b 1 (t) dT b dt + b 2 (t), (23)
the time functions b 1 (t) and b 2 (t) are adapted to give a second-order ordinary differential equation with constant coefficients for T b -T t b whose response is set to the critically damped regime. It is then found that

b 1 (t) = - 2ρ b c p b τ b 2 (t) = ρ b c p b T t b -T b (t) τ 2 + 1 2δ d dt q w,1 +q w,2 ρ b c p b -d dt 1 ρ b c p b V V Qtot dV (24)
Different choices are possible for the initial condition Ω(t = 0). It is set to a reference value Ω re f = q re f w δ which is estimated from empirical Nusselt number formula for the desired operating condition.

Validation

Temperature as a passive scalar

Temperature is first treated as a passive scalar in order to validate the proposed method in the case of a flow with constant thermo-physical properties. The bulk Reynolds number is then not sensitive to the variations of temperature and is set as Re b = 5850. The wall temperature are identical T w,c = T w,h = 1000 K and the initial instantaneous field of temperature is obtained for a channel flow simulation where T b = 1100 K. The target bulk temperature is set to T t b = 2000 K. The wall heat flux averaged in homogeneous directions on both walls can be expressed as a Nusselt number Nu defined as

Nu = 4q w δ λ w (T b -T w ) , (25)
where the reference length scale 4δ corresponds to the channel flow hydraulic diameter and where λ w is the thermal conductivity at the walls. Figure 4 (a) shows for two different decay times τ that the bulk temperature follows the theoretical response curve. Similarly, to the wall shear stress, for small imposed value of τ, the average wall heat flux shown in figure 4 (b) needs an incompressible physical time to adapt to the change of T b . As expected, the obtained steady value of the wall heat flux is larger than its initial value due to the significant change in the bulk temperature. The Nusselt number is shown in figure 4 (c). Its final value remains close to the initial one because of the weak dependency of Nu with the parameter T b /T w .

Simultaneous control of bulk Reynolds number and bulk temperature

The temperature field is now coupled to the velocity field by taking into account the change of properties with temperature. Applying the control method for both Re b and T b is necessary as controlling and varying one quantity τ . The temporal evolution of bulk and friction Reynolds numbers are presented in figure 5. Deviation of the bulk Reynolds number is larger compared to results in section 2.2.2 due to the variation of the bulk temperature from 1100 K to 2000 K which induces large temporal changes in physical properties such as the bulk dynamic viscosity. Despite such strong perturbations, the controlling approach derived for constant properties flow remains robust and brings Re b to its target within a duration close to the desired decay time.

Similarly, the bulk temperature shown in figure 6 (a) deviates slightly from the theoretical response curve before reaching a plateau at the desired target value T t b =2000 K. Corresponding wall heat flux and Nusselt number are given in figure 6 (b) and (c), respectively. This time, the final value taken by the Nusselt number changes because of the modified bulk Reynolds number.

This coupled case with simultaneous control of the bulk Reynolds number and the bulk temperature demonstrates the adequacy of the method even with strong deviations from the assumptions made to derive the controllers. As mentioned previously, the equations to derive the method can be modified to account for the temporal evolutions of bulk variable properties ( see next section). As shown here, relying instead on the robustness of the obtained constant-properties controllers is also a good strategy.

Conclusion

Bulk Reynolds number and bulk temperature are key quantities when reporting results in channel flow simulations. There are situations are when one wishes to control these parameters although they are actually results of the simulations when imposing a fixed source term in the momentum and energy equations. A control strategy has been developed to carefully set the desired target values for the bulk Reynolds number and bulk temperature by temporally adapting the prescribed source terms. Final equations appear as PI controllers with a thorough setting of the coefficients so that the dynamics of the quantities of interest is exactly controlled when the flow thermo-physical properties are constant or even variable. The controlling method is validated in several cases and is shown to be robust and efficient. Such a method enables to safely study the effect fo different numerical or physical description of the flow while safely maintaining identical bulk Reynolds numbers and/or bulk temperatures.

Figure 1 :

 1 Figure 1: Channel flow configuration of half-width δ. In heat transfer simulations, both wall temperatures T w,c (lower wall) and T w,h (upper wall) are prescribed where T w,h ≥ T w,c .

Figure 2 :

 2 Figure 2: Temporal evolution of Re b (a) and Re τ (b) for an isothermal case (T =1000K). Circles: theoretical solution; Black plain line: τ = 0.33δ/u re f τ ; Gray plain line: τ = 0.033δ/u re f τ ; Black dashed line: τ = 0.0033δ/u re f τ .

Figure 3 :

 3 Figure 3: Temporal evolution of Re b (a), Re τ on the cold and hot sides (b) for a non-isothermal case (T w,c =950K,T w,h =2050K). Circles: theoretical solution; Black plain line: τ = 0.33δ/u re f τ ; Gray plain line: τ = 0.033δ/u re f τ .

Figure 4 :

 4 Figure 4: Temporal evolution of T b in Kelvin (a), q w in W/m 2 (b) and Nu for Re b = 5850 with contant properties. Circles: theoretical solution; Black plain line: τ = 0.33δ/u re f τ ; Gray plain line: τ = 0.033δ/u re f τ .

Figure 5 :Figure 6 :

 56 Figure 5: Temporal evolution of Re b (a) and Re τ (b) for with varying T b at the same time. Circles: theoretical solution; Black plain line: τ = 0.33δ/u re f τ ; Gray plain line: τ = 0.033δ/u re f τ .

Figure 8 :

 8 Figure 8: Temporal evolution of T b in Kelvin (a) and Nusselt number (b) for varying Re b at the same time with the control equations. Circles: theoretical solution; Black plain line: τ = 0.33δ/u re f τ ; Gray plain line: τ = 0.033δ/u re f τ .

Acknowledgements

This work was granted access to the HPC resources of [CCRT/CINES/IDRIS] under the allocations 2013-2b0164 and 2015-2b0164 made by GENCI (Grand Equipement National de Calcul Intensif). Vincent Moureau and Ghislain Lartigue from the CORIA lab, and the SUCCESS scientific group are acknowledged for providing the YALES2 code.

Accounting for variable flow properties during control of bulk Reynolds number and bulk temperature

For the sake of completeness, controllers that account for temporal variations of thermo-physical properties of the mixture are derived. The equations 7 and 21 and for the derivative of Re b and T b are then replaced by

Computing the temporal derivative of theses expressions gives modified expressions for the control equations. Hence, the momentum source term is still set by the differential equation dS 1 dt = a 1 (t) dRe b dt +a 2 (t) but with modified coefficients:

Similarly, it is found that

with the energy source term determined by dΩ dt = b 1 (t) dT b dt + b 2 (t). The same simulation case as in section 3.2.2 is carried out with these new equations. The obtained temporal evolutions of bulk Reynolds number and temperature are shown in figures 7 and 8 with their corresponding profiles of Re τ and Nusselt number. Thanks to the modified equations to take into account temporal variations of thermo-physical properties, the bulk quantities henceforth follow their prescribed behavior.


Sylvain Leclercq
email: sylvain.leclercq@edf.fr

Gilles Rousselier

G Cailletaud

A generic method for modeling the behavior of anisotropic metallic materials : application to recrystallized zirconium alloys

Keywords: anisotropic material, constitutive behaviour, crystal plasticity, polycrystalline material, Zirconium alloys

come

A generic method for modelling the behavior of anisotropic metallic materials : application to recrystallized zirconium alloys

Introduction

Since the pioneering work of [START_REF] Hill | The Mathematical Theory of Plasticity[END_REF], several efforts have been made to improve the description of initial plastic anisotropy according to macroscopic constitutive equations (see for instance, [START_REF] Gotoh | A theory of plastic anisotropy based on a yield function of fourth order (plane stress state) -I/II[END_REF][START_REF] Barlat | A six-component yield function for anisotropic materials[END_REF][START_REF] Karafillis | A general anisotropic yield criterion using bounds and a transformation weighting tensor[END_REF][START_REF] Bron | A yield function for anisotropic materials -Application to aluminum alloys[END_REF]. Nevertheless, these macroscopic models generally fail to accurately represent non proportional and cyclic loadings [START_REF] Cailletaud | Utilisation de modèles polycristallins pour le calcul par éléments finis[END_REF]. Especially, the simulation and prediction of the yield surface distorsion remain difficult to achieve with classical phenomenological models, and more or less complicated improvements of these models do not give total satisfaction [START_REF] Vincent | An improvement of multiaxial ratchetting modelling via yield surface distortion[END_REF][START_REF] Kuroda | A phenomenological plasticity model with nonnormality effects representing observations in crystal plasticity[END_REF].

Polycrystalline models generally give a good description of the plastic anisotropy and yield surface distorsion, because they are based intrinsically on slip lattice as well as texture informations [START_REF] Cailletaud | A micromechanical approach to inelastic behavior of metals[END_REF][START_REF] Pilvin | Approches multiéchelles pour la prévision du comportement anélastique des métaux[END_REF][START_REF] Calloch | Essais triaxiaux non-proportionnels et ingéniérie des modèles de comportement[END_REF]. Nevertheless, their main disadvantage is the large computation times needed for parameter identification and finite element calculations. Indeed, these models are built in such a way that the number of internal variables to compute may be a thousand times larger than the number of internal variables using a phenomenological model. Thus one understands the relatively rare use of polycrystalline models in industrial calculations.

In the field of nuclear industry, important work has been done on the viscoplastic behavior of zirconium alloys, using phenomenological models as well as polycrystalline approaches, both combined to experimental multiaxial tests performed on tubes. These multiaxial tests (pure tension, biaxial tension, tensiontorsion) are necessary to investigate the material response to the complex loadings that are applied to the nuclear fuel rods in Pressurized Water Reactors. It can be seen that an accurate description of the anisotropy and irradiation effects requires quite a complex set of equations in phenomenological approaches [START_REF] Robinet | Etude expérimentale et modélisation du comportement viscoplastique anisotrope du Zircaloly 4 dans deux états métallurgiques[END_REF][START_REF] Schaeffler | Modélisation du comportement élasto-viscoplastique anisotrope des tubes de gaine du crayon combustible entre zero et quatre cycles de fonctionnement en réacteur à eau pressurisée[END_REF][START_REF] Richard | Modeling of the cold work stress relieved Zircaloy-4 cladding tubes mechanical behavior under PWR operating conditions[END_REF] while the multi-scale approach allows much simpler developments [START_REF] Lebensohn | A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycristals: application to zirconium alloys[END_REF][START_REF] Geyer | Etude et modélisation du comportement élastoviscoplastique de tubes en Zircaloy par une approche polycristalline[END_REF][START_REF] Onimus | Approche expérimentale et modélisation micromécanique du comportement des alliages de zirconium irradiés[END_REF]. Nevertheless, none of these models is used at present time in really industrial finite element calculations because of their complexity or too large CPU time.

In order to combine the good predictive ability of polycrystalline approaches and the moderate CPU time necessary to conduct safety analyses for the nuclear industry, we propose in the present paper a simplified polycrystalline model (the so-called RL model - Rousselier and Leclercq, 2004) to simulate the anisotropic viscoplastic behavior of metallic materials. A generic method is presented that makes it possible to build a simplified anisotropic material texture, based on the principal features of the pole figures. The method is applied to recrystallised zirconium alloys. The model aims at representing the mechanical behavior of polycrystals with a minimum number (between 5 and 10) of crystalline orientations ("grains"), and considers only 6 slip systems by grain for describing the inelastic strain tensor, whatever the crystallographic lattice of the material under study is. It must be very clear that our goal is not to represent accurately the material texture, but only its mechanical effect on the response of a structure to a mechanical loading. Nevertheless, we show in the following that our choice of the crystallographic orientations is closely related to the real material texture, and moreover that it can be identified thanks to an optimization algorithm by comparison with the results of experimental tests.

In the next section, we recall the major features of the simplified viscoplastic RL-polycrystalline model. In section 3, we show how the model has been applied to the simulation of the anisotropic behavior of a recrystallized zirconium alloy. Particular attention is paid to the description of the parameter identification of the model using several multiaxial tests. An example of finite element calculation performed on a tube is presented in section 4. Concluding remarks are given in the last section of this paper.

The simplified viscoplastic RL-polycrystalline model

The model is presented in the case of small deformations and small rotations. We want to emphasize that the RL model does not depart from classical polycrystalline models. Consequently, all existing and future theoretical and software developments apply easily to the model, with regard for example to finite deformations, rotation of crystallographic orientations, localization self-consistent models, or slip systems constitutive equations.

Let s n be the normal vector of a slip plane and s l the unit slip direction vector. The orientation tensor of the slip system number s is defined as:

2 / ) ( s s s s s n l l n m ⊗ + ⊗ =
(1)

We consider only 6 universal slip systems corresponding to shear and extension strain rates in the local frame of a given crystallographic orientation (Figure 1):

3 2 1 x x x [ ] 0 0 1 1 = n , [ ] 0 1 0 2 = n ,
[ ]

1 0 0 3 = n [ ] 0 1 0 1 = l ,
[ ]

1 0 0 2 = l , [ ] 0 0 1 3 = l ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = 0 0 0 0 0 1 0 1 0 2 1 1 m , ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = 0 1 0 1 0 0 0 0 0 2 1 2 m , ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = 0 0 1 0 0 0 1 0 0 2 1 3 m [ ] 0 2 / 1 2 / 1 4 = n , [ ] 2 / 1 2 / 1 0 5 = n , [ ] 2 / 1 0 2 / 1 6 = n [ ] 0 2 / 1 2 / 1 4 - = l , [ ] 2 / 1 2 / 1 0 5 - = l , [ ] 2 / 1 0 2 / 1 6 - = l ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡- = 0 0 0 0 1 0 0 0 1 2 1 4 m , ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - = 1 0 0 0 1 0 0 0 0 2 1 5 m , ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ - = 1 0 0 0 0 0 0 0 1 2 1 6 m
(2)

The incompressible plastic strain rate of a given orientation or "grain" results from the multiple slip:

⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - - - = = ∑ = 6 5 2 3 2 5 4 1 3 1 4 6 6 1 2 1 γ γ γ γ γ γ γ γ γ γ γ γ γ ε
For hcp crystals, the deformation results from 5 main slip families (for instance, see [START_REF] Geyer | Etude et modélisation du comportement élastoviscoplastique de tubes en Zircaloy par une approche polycristalline[END_REF][START_REF] Starolesky | A constitutive model for hcp materials deforming by slip and twinning: application to magnesium alloy AZ31[END_REF]: basal B<a>, prismatic P<a>, pyramidal π1<a>, π1<c+a> and π2<c+a> (Figure 2a). Only pyramidal slips <c+a> contribute to the strain rate component . According to equation (3), the systems 5 and 6 of the model will be associated with these physical slip systems.

For the same reason, the systems 2 and 3 of the model will be associated with basal slips that contribute to the strain rate components and . Finally, the systems 1 and 4 of the model will be associated with prismatic slips that contribute to the strain rate components , , and . So, for hcp materials, we should consider 3 slip families (1,4), (2,3) and (5,6) that will have the properties respectively of prismatic, basal and pyramidal slips (Figure 2b). In order to take into account the fact that systems 1 and 4 have not the same schematic representation in Figure 2b, we choose to split the "prismatic family" into two parts, and we assume that they have not the same properties, especially in terms of critical resolved shear stress. We do not pretend that these non physical slip systems accurately represent the physical ones. Again it must be clear that our goal is only to model the macroscopic mechanical behavior of the material. Any constitutive equation for the slip systems can be used with the RLpolycrystal model, for example the power law form of [START_REF] Asaro | Texture development and strain hardening in rate dependent polycrystal[END_REF] or equations depending on dislocation densities (for instance, see [START_REF] Hoc | Polycrystal modeling of IF-Ti steel under complex loading path[END_REF][START_REF] Hoc | A procedure for identifying the plastic behavior of single crystals from the local response of polycrystals[END_REF]. In the following, we will use a phenomenological model with two internal variables, isotropic and kinematic , as proposed by [START_REF] Cailletaud | A micromechanical approach to inelastic behavior of metals[END_REF] :

s r s X ( s s s s X Sign v -) = τ γ & & (4) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - = s n s s s s s K r X Max v τ , 0 & (5) ( ) [ ] [ s t t t st c s s bv Q v b Q r - - + - - + = ∑ = exp( 1 exp 1 6 1 τ ] (6) s s s s s s v X c a X & & & - = γ (7)
Coefficients are the components of the hardening moduli matrix st Q Q . The diagonal terms represent the self-hardening of each system, and the non-diagonal terms the latent hardening. Generally speaking, Q can be non-symmetric. For simplicity, in the following, we have chosen a symmetric form for Q . This induces 8 different values to take into account the 4 slip families of the hcp structure. There is no latent hardening for the kinematic variables. is the initial critical value of the resolved shear stress c s τ s τ on the slip system s, given by:

s g s m : σ τ = (8)
where g σ is the stress tensor in the "grain" number g (from 1 to ). This stress tensor N g σ is evaluated using localization equations. In the same way as for the constitutive equations of the slip systems, several choices can be made. One can refer to the general presentation of the RL model (Rousselier and Leclercq, 2004) to have an overview of what may be used in the framework of a simplified polycrystalline approach. In the present application, we consider the localization equation proposed by [START_REF] Cailletaud | Une approche micromécanique phénoménologique du comportement inélastique des métaux[END_REF][START_REF] Cailletaud | A micromechanical approach to inelastic behavior of metals[END_REF], in which auxiliary strain variables g β are memorized for each "grain" :

) ( g g B β α σ - + Σ = with ∑ = = N g g g f B 1 β (9) ) 1 ( 2 β μ α - = and ν ν β - - = 1 5 4 15 2
where μ and ν are the shear modulus and the Poisson coefficient, respectively. In equation ( 9), Σ is the macroscopic stress tensor (over the Representative Volume of Element) and the volume fraction of "grain" g. g f

The evolution of g β is given by a phenomenological kinematic equation, with 2 adjustable parameters [START_REF] Cailletaud | Une approche micromécanique phénoménologique du comportement inélastique des métaux[END_REF] and [START_REF] Pilvin | Approches multiéchelles pour la prévision du comportement anélastique des métaux[END_REF]:

D * D ( ) 3 2 2 * p g p g g p g g J D D ε ε β ε β & & & ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - - = (10)
where

( ) ( ) 2 / 1 2 2 / : d d T T T J
= and d T is the deviator of a tensor T .

In the present study, the parameter is chosen to be zero.

*

D

Let us note the phenomenological character of the present scale transition rule. The idea behind the formulation is that the development of a plastic accommodation will decrease the level of the residual intergranular stresses. The shape of the expression is imported from the classical macroscopic models. The original authors (Pilvin and Cailletaud, 1990) suggest to calibrate the parameters of the rule on finite element computations of a representative material element.

This has not been made in the present study where, following the spirit of the model proposed in this paper, we have decided to consider the variables in the set of parameters to be identified. Nevertheless, the values found in the direct identification process are in good agreement with classical values of these parameters.

Application to the simulation of the recrystallized Zircaloy behavior

Material

The material under study consists in Zircaloy-4 tubes, the weight chemical composition of which is reported in Table 1.

The material is tested in the α-recrystallized state, i.e. after processing, the tubes are submitted to a heat treatment of 700°C ± 30°C during 4 to 5 hours. After recrystallization, the microstructure consists in small equiaxed grains (2 to 5 μm)

which are bounded with larger ones (15 to 18 μm) (Bouffioux, 2004).

The reference frame of the material is the frame

3 2 1 X X X z rθ of the tubes.
The texture of these tubes is strongly anisotropic, see Figure 3 [START_REF] Robinet | Etude expérimentale et modélisation du comportement viscoplastique anisotrope du Zircaloly 4 dans deux états métallurgiques[END_REF]. It Note that the texture presented here results in the combination of two effects : a material one (due to the hcp lattice, and to the heat treatment), and an effect due to tube processing. We have to point out that is texture is characteristic of tubes made of recrystallized zirconium alloys (tubes made of cold-word stress relieved material would have significantly different pole figures).

Experimental data base and parameter identification

Several experimental tests have been performed by [START_REF] Robinet | Etude expérimentale et modélisation du comportement viscoplastique anisotrope du Zircaloly 4 dans deux états métallurgiques[END_REF]

- - = s E & 1 7 10 . 6 , 6 - - = s E & 1 4 2 10 . 6 , 6 ) ( - - = = s E J E & &
). The validation data base (i.e. the data base that is used to demonstrate the predictive ability of the model) is built with the remaining experimental tests.

Usually three types of parameters are identified when one uses a polycrystalline approach :

-the material parameters linked with the slip properties of each system (critical resolved shear stress, viscosity, hardening,…) ;

-the matrix defining self-and latent-hardening ( ) ; We have considered α to be a parameter, instead of using its theoretical value ) 1 ( 2

β μ α - =
, Eq. ( 9) .

In the usual case, the texture of the material is an input data which is made of a certain number of grains represented by their crystallographic orientation.

Generally, the number of grains is large. For example, [START_REF] Geyer | Etude et modélisation du comportement élastoviscoplastique de tubes en Zircaloy par une approche polycristalline[END_REF] Identifying these parameters gives us the opportunity to verify the ability of the model to match not only with macroscopic tests, but also with the real material texture.

Thus, the parameters of the RL-polycrystal model have been identified using the identification data base defined above. The parameter values are reported in Table 3. As mentioned above, four slip families are considered (i.e., 1, 4, 2-3 and 5-6), each family having its own set of parameters. The self-and latent-hardening matrix has the following form : 2). A second validation of this original choice will follow with the study of the predictive ability of the RL model applied to the recrystallized Zircaloy-4.

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = 5
Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q ( 11

Predictive ability of the model and discussion

In order to investigate its predictive ability, the model has been tested on several loading situations which had not been used in the identification data base.

These multiaxial tests are conducted at the same equivalent (in the von Mises sense) strain rate (

1 4 10 . 6 , 6 - - = s E &
). We present here three tension-shear tests with three different strain biaxiality ratios (

zz z E E θ
), and two tension-internal pressure tests, with two different strain biaxiality ratios (

zz E E θθ
). The simulated results

are compared with the experimental ones on Figures 7 to 11.

Several features are to be pointed out:

i. First of all, we have to precise that our goal was to test the RL model ii.

As concerns the RL model, agreements between simulation and experiments are good, for several loading situations and especially the tension-internal pressure tests which were not represented in the identification data base. This confirms the ability of the model to represent situations which have not been taken into account in the identification procedure, especially as concerns anisotropy. We would like to point out that getting this result was not obvious, because of the small number of grains used to simulate the texture of the Zircaloy-4 tube.

iii.

An important feature, linked with the previous one, is that the RL model can be totally identified with only two experimental tests, representing two loading directions (pure tension , and pure shear zz z θ ). A confirmation of this can be found on Figures 12 and13, where the critical stresses of the 6 slip systems are represented for each grain (from grain 1 to grain 6), in the macroscopic -zz θθ and -zz θ z planes. It can be seen here that the use of a tension test conducted up to minimum 500 MPa allows to activate and thus identify the material parameters of slip systems # 1, 4, 5 and 6. Note that in the case of a tensile test, slip systems # 2 and 3 are never activated in the RL model. It is interesting to remember that these two systems are associated to basal slip, which is known to be rarely observed when studying the experimental mechanisms of unirradiated Zircaloy-4 [START_REF] Geyer | Etude et modélisation du comportement élastoviscoplastique de tubes en Zircaloy par une approche polycristalline[END_REF]. Using a pure shear test conducted up to 250 MPa allows one to activate slip systems # 1, 4 and 3. As the material parameters of slip system # 2 are equal to those of slip system # 3, it is therefore demonstrated that every slip system has been activated at least once in the simulation of the identification data base.

This explains why these two testing directions are sufficient for the parameter identification, provided the experiments are conducted far enough in terms of strain (and thus stress) value. The activation of many slip systems at ≈500 MPa in Figure 12 explains why the pressure curves are so flat in Figures 10 and 11, a feature that is not due to curve fitting but to the texture of the material and cannot be predicted by macroscopic models. v.

The present study shows that, in the case of an anisotropic material such as Zircaloy, the parameters of the RL model can be identified knowing the general features of the material texture, and with experimental tests performed in two directions (pure tension and pure shear) provided these tests are conducted sufficiently far so as to activate each slip system. With these informations, it is possible to get material parameters that are robust enough to ensure a good predictive ability, which is not the case for macroscopic approaches. Of course, enlarging the identification data base with tests performed in other directions allows to improve the parameter identification of both kind of models, but at a higher cost. Once again, introduction of other experimental tests would be essential to identify parameters linked with specific material features such as creep or cyclic effects, which are not the concern of this paper.

Another way to verify the robustness of the material parameters (and of the proposed constitutive equations) is to investigate whether they are dependent on the identification data base or not. In order to answer this question, we have added to the identification data base the equibiaxial tension with internal pressure test

( 1 = zz E E θθ
). With the use of this extended data base, we have identified the material parameters of the RL model from the same initial values as with the first identification. The results are presented in able 4. It is remarkable to point out how the values of these new parameters are close to those presented in Table 3. This feature is another proof that a relatively poor identification data base is sufficient to get robust material parameters of the RL model, and that this model is intrinsically predictive.

A finite element calculation

As the principal purpose of this paper is to propose a material model based on the classical polycrystalline approach and compatible with finite element calculations (i.e., consuming "reasonable" CPU time), we present in the following the simulation of a Zircaloy-4 tube submitted to an internal pressure and an imposed axial displacement. The structure is calculated using both RL model and the Chaboche+Hill macroscopic approach. We compare the numerical results in terms of strain-stress curves, as well as the CPU time needed to perform the calculation.

The 3D mesh of the finite element calculation is presented Figure 14. It consists in a ring, part of a Zircaloy-4 tube, meshed with 200 quadratic elements with 20 nodes. The imposed boundary conditions are such that the axial strain of the ring is zero, and the hoop strain rate is equal to 6,6 10 -4 s -1 . Two finite element calculations have been performed, with the Chaboche+Hill approach and the RL model. The macroscopic strain and stress states in the inner part of the tube are reported Figure 15, compared with the experimental data. The results are consistent with the ones shown on Figure 11, which shows the 0D (representative elementary volume) simulation of this test. It can be pointed out that, especially for the tensile results in Figures 11 and15, the smaller strain behavior is better predicted with the Chaboche+Hill approach, and the RL model fits better for the larger strains. One explanation for this may be found in the fact that the Chaboche+Hill approach is written in a purely viscoplastic way (the viscoplastic strain always exists, even if it is very small for small total strains). On the contrary, the RL model is written in such a way that viscoplastic slip is activated on a given system if and only if the critical value of the resolved shear stress is reached. Thus non elastic strain is intrinsically smaller in the RL approach than in the classical one. At higher strains, the problem is different : in the Chaboche+Hill model the 7

(one tensorial kinematic and one scalar) hardening internal variables saturate rapidly (relatively low strain hardening) and in the RL model the possible activation of 36 (6 grains x 6 slip systems) slip systems makes it possible to have a better description of the global material behavior.

For complete understanding, we have to note that the hoop stress gradient between the inner and the out part of the tube (tube thickness = 570 micrometers, inner radius = 480 millimeters) is approximately 80 MPa. This confirms that the tube may be considered as a real structure, and not as a thin tube.

Concluding remarks

In the present paper, a simplified polycrystalline model (the so-called RL model) is proposed to simulate the anisotropic viscoplastic behavior of metallic materials. A generic method is presented that makes it possible to build a simplified anisotropic material texture, based on the principal features of the pole figures. The method is applied to a recrystallised zirconium alloy, used as clad material in the fuel rods of nuclear power plants. The objective of this study was to propose a robust, predictive but little CPU time consuming model for industrial finite element calculations.

In order to achieve this goal, several simplifications of the classical polycrystalline approach have been performed : the definition of slip systems and the choice and number of crystalline orientations are the major aspects developed here. As a matter of fact, these simplifications lower the physical basis of the approach proposed in this paper, but it has been shown that the main characteristics of the classical polycrystalline approach were kept : predictive ability and intrinsic way to describe the material anisotropy. Moreover, we have shown that the simplicity of the RL model allows one to follow a parameter identification scheme that can be dedicated to each slip system. The relatively small CPU time (as compared to standard polycrystalline approaches), combined with the possibility to follow each slip system, makes it possible to get robust material parameters with a minimum number of tests in the identification data base.

The finite element calculations performed with the RL model confirm the promising features of this model in terms of good quality of the simulations, and CPU time of the same order as with more classical macroscopic models.

The next step, under progress, is to apply the material parameters of the model to another texture, obtained from recrystallized zirconium alloy sheets, of the same chemical composition. The results, combined with the modeling of zirconium alloys' damage properties, will be exposed in a forthcoming paper. ).

Experimental points, RL model (solid line), Chaboche+Hill (bold line). ).

Experimental points, RL model (solid line), Chaboche+Hill (bold line). ).

Figure 15. Macroscopic strain and stress states at the outer part of the ring.

Experimental points, RL model (solid line), Chaboche+Hill (bold line).

  shows z θ pole figure { 0002 }, i.e. basal poles corresponding to axis c of the crystal or axis , and pole figure i.e. prismatic poles corresponding to axis in RL model. 1 x

  st Q -the parameters of the localization equation ( and D α ).

Figure

  Figure 4a presents the pole figures (axis x3-left and axis x1-right) of the six

)

  One has to point out that, even if the RL slip families are not physical ones, the identified values of the critical resolved shear stresses are in agreement with what is classically observed on non-irradiated zirconium alloys, i.e.c τ prismatic < c τ basal < c τ pyramidalThe numerical values obtained for the whole set of parameters are consistent with the ones obtained with standard polycrystalline approaches on the same type of material (see for instance[START_REF] Geyer | Etude et modélisation du comportement élastoviscoplastique de tubes en Zircaloy par une approche polycristalline[END_REF].The comparisons between experiments and simulations of the identification data base are presented on Figs.5 and 6. The simple lines present RL model simulations. The bold lines correspond to simulations performed with the use of a classical macroscopic viscoplastic model (using one kinematic and one isotropic hardening, see for instance[START_REF] Lemaitre | Mechanics of solid materials[END_REF] in which the yield surface is defined by a Hill criterion, in order to take into account the Zircaloy anisotropy. Of course, this model has been identified with the same data base as the RL model. Use of this model is made in order to have a comparison basis between the approach proposed in this paper and more classical models still widely used. The simulated results are good, and cannot allow one to distinguish between the two approaches.On Figure4b, one can see the pole figures (axis x3-left and axis x1-right) of the 6 grains representing the texture of the material under study, after identification of the 3 constants a ψ , b ϕ , c ϕ . These pole figures are to be compared with the ones presented on Figure 3. The good agreement observed here is a first validation of our choice to integrate 3 Euler angles into the parameter identification procedure (see Table

  as concerns its ability to represent the anisotropic features of the material in terms of yield stress and flow direction. To this aim, we have used only constant strain rate tests to identify the model parameters. As a matter of fact we do not, in the following, draw any conclusion on the predictive ability of the response of the model to other mechanical loadings cyclic tests, creep or relaxation tests, … Particularly, we are aware that only two experimental tests are not sufficient enough to identify the model parameters associated with all the material characteristics.

  the macroscopic approach shows the classical feature that the Hill yield surface cannot be identified without appropriate experimental tests in the three directions. Typically, our identification data base does not include any internal pressure test, and this results in the impossibility of identifying any anisotropy between the rr , θθ , and directions for the macroscopic model, as shown on Figs. 10 and 11. As a matter of fact, the identification data base makes it possible to identify only the shear ( zz z θ ) anisotropy for this kind of model.

Figure captions Figure 1 .

 captions1 Figure captions

Figure 2 .

 2 Figure 2. 2a : the physical slip families of a hcp crystal. 2b : the 4 slip

Figure 3 .

 3 Figure 3. The experimental texture of the recrystallized Zircaloy-4. 3a : the

Figure 4 .

 4 Figure 4. The pole figures of the simulated material (axis x3-left and axis

Figure 5 .

 5 Figure 5. Identification data base : tension tests. Experimental points, RL

Figure 6 .

 6 Figure 6. Identification data base : pure shear test. Experimental points, RL

Figure 7 .

 7 Figure 7. Simulation data base : tension-shear ( 4 3 = zz z E E θ).

Figure

  Figure 8. Simulation data base : tension-shear ( 2 3 = zz z E E θ).

Figure

  Figure 9. Simulation data base : tension-shear ( 3 = zz z E E θ). Experimental

Figure 10 .

 10 Figure 10. Simulation data base : tension-internal pressure ( 1 = zz E E θθ).

Figure 11 .

 11 Figure 11. Simulation data base : tension-internal pressure (

Figure 12 .

 12 Figure 12. Critical stresses of the 6 slip systems represented for the material

Figure 13 .

 13 Figure 13. Critical stresses of the 6 slip systems represented for the material

Figure 14 .

 14 Figure 14. 3D mesh of a ring submitted to tension and internal pressure ( ∞ = zz E E θθ).

  Figure 2

  Figure 4b

  Figure 14

p 33 ε& p 13 ε& p 23 ε& p 11 ε& p 12 ε& p 22 ε&

  Table V presents the CPU times of the two finite element calculations. One can see that CPU time of RL model is not much larger than the one of Chaboche+Hill model (ratio of approximately 5), compared with the accuracy of the simulated results. If one recalls that the Chaboche+Hill model used in the present study needs to integrate 13 internal variables (the elasticity tensor, one tensorial kinematic hardening, one scalar isotropic hardening, the equivalent viscoplastic strain), and that the RL model requires the integration of 48 internal variables, the ratio of 5 in the CPU times obtained in Table5can be easily related to the ratio of the number of internal variables to integrate. It must also be pointed

	out that the same finite element calculation performed with a standard
	polycrystalline model (for instance Geyer, 1999, with 7206 internal variables to
	integrate) would have last approximately 200 times longer than the one with RL
	model.

Table 1 .

 1 Weight chemical composition (%) of the Zircaloy-4 under study.

	Cr	Fe	Sn	O	Zr
	0,10	0,21	1,3	0,135	balance
	Al (ppm)	C (ppm)	H (ppm)	Hf (ppm)	N (ppm)	Si (ppm)
	98	151	9	56	24	99

Table 2 .

 2 Euler angles of the 6 grains representing the recrystallized Zircaloy-4 texture.

Table 3 .

 3 Material parameters of the RL model identified on the recrystallized Zircaloy-4, with the initial identification data base.

	ψ	a	ϕ	b	ϕ	c	α (MPa) D
	54° E (MPa) ν 65°	73,7°	1 Q (MPa)	37963 2 Q (MPa)	256 5 Q (MPa)
	125000	0.36				77,4	88,8	78
	12 Q (MPa)	14 Q (MPa) 15 Q (MPa)	23 Q (MPa)	56 Q (MPa)
	198,7	116,3	111	47,7	175
	Q (MPa)	b
	9,8	1,14
	c 1 τ (MPa)	1 n	1 K (MPa)	c 2 τ (MPa)	n	2	2 K (MPa)
	31,4	1,34	8730	227,6	1,63	5036
	c 4 τ (MPa)	n	4	4 K (MPa)	c 5 τ (MPa)	5 n	5 K (MPa)
	146,2	1,12	1396,3	231,7	1,07	3463,5
	1 a (MPa)	1 b	1 c	2 a (MPa)	2 b	c	2
	5400,16	5,81	920,89	2358,34	7,77	641,28
	4 a (MPa)	4 b	c	4	5 a (MPa)	5 b	c	5
	5321,46	6,33	765,4	6522,16	2,74	972,26

Table 4 .

 4 Material parameters of the RL model identified on the recrystallized Zircaloy-4, with the extended identification data base.

	ψ	a	ϕ	b	ϕ	c	α (MPa) D
	57° E (MPa) ν 66,7°	67,3°	1 Q (MPa)	37633 2 Q (MPa)	260 5 Q (MPa)
	125000	0,36				78,4	89,5	76,4
	12 Q (MPa)	14 Q (MPa) 15 Q (MPa)	23 Q (MPa)	56 Q (MPa)
	195,8	116,3	112	49	171,3
	Q (MPa)	b
	9,4	1,2
	c 1 τ (MPa)	1 n	1 K (MPa)	c 2 τ (MPa)	n	2	2 K (MPa)
	32,4	1,37	8879	231,7	1,65	5033
	c 4 τ (MPa)	n	4	4 K (MPa)	c 5 τ (MPa)	5 n	5 K (MPa)
	141,2	1,12	1394,1	242,1	1,07	3471,5
	1 a (MPa)	1 b	1 c	2 a (MPa)	2 b	c	2
	5415,1	5,7	908,88	2343,29	7,8	628,53
	4 a (MPa)	4 b	c	4	5 a (MPa)	5 b	c	5
	5262,83	6,24	765,9	6641,55	2,69	972,26

Table 5 .

 5 CPU times (PC/Linux) obtained with Chaboche+Hill and RL models.

	Model	Hill	RL
	CPU time (s)	196	977

& & & & & & & & & & & & & &

Acknowledgements

The authors would like to thank Dr. Jacques Besson (Centre des Matériaux, Ecole des Mines de Paris) for his fruitful help in the finite element calculations.