H Q Nguyen

L Gallimard

C Bathias

Numerical simulation of the coupling between thermal dissipation and fish-eye crack growth in very high cycle fatigue regime

Keywords: fatigue crack growth, finite element model, fish-eye fracture, plasticity, thermomechanical analysis, ultrasonic fatigue, very high cycle fatigue

A C T In this paper, we study the temperature field associated with the propagation of a fatigue crack in a very high cycle fatigue regime during ultrasonic fatigue testing. We use a Paris-Hertzberg crack growth law to compute the evolution of the crack and a perfectly elastic-plastic constitutive law to compute the plastic dissipation per cycle at the tip of the crack. A thermomechanical finite element model is proposed to estimate the evolution of the temperature field during the crack propagation. Numerical results obtained agree fairly well with experimental results.

The considerable attention given, in the past ten years, to fatigue fracture of metallic material in very high cycle regime, is mainly due to the increasing use of engineering materials in applications with service lives reach-Correspondence: L. Gallimard. E-mail: laurent.gallimard@u-paris10.fr ing up to 10 10 load cycles. It has been shown that very high cycles fatigue (VHCF) failures of steels are characterized by a fracture which does not occurs on the surface but rather internally in the material, especially in high strength steels. [START_REF] Bathias | How and why the fatigue sn curve does not approach a horizontal asymptote[END_REF][START_REF] Sohar | Gigacycle fatigue behavior of a high chromium alloyed cold work tool steel[END_REF][START_REF] Li | Effect of loading type on fatigue properties of high strength bearing steel in very high cycle regime[END_REF] The initiation of fatigue cracks can be considered differently from physical or mechanical points of view. At the microscopic level, Mughrabi shows that the initiation of fatigue crack in the gigacycle fatigue regime can be described in terms of a micro structurally irreversible portion of the cumulative cycle strain. This means that there is no basic difference between fatigue mechanisms in low, mega, and gigacycle fatigue except for the strain localization. However, specific mechanisms can occur depending on the fatigue life. The fatigue life seems to be a key parameter to determine correctly the fatigue initiation location. In low cycle fatigue, in mega cycle fatigue and in gigacycle fatigue different mechanisms can operate at different scales of plasticity. In the low cycle fatigue regime, the cyclic plastic deformation is critical at the surface but exists also in the bulk of the metal. Typically, several cracks nucleate from the surface. When the fatigue life is below 10 5 cycles, general plastic deformation of the specimen bulk governs the initiation. When the fatigue life is between 10 6 and 10 7 cycles the plastic deformation depends on the plane stress surface effect and the presence of flaws, which explain the critical location of fatigue initiation. Typically, the initiation starts with one crack only, from the surface. However approaching 10 9 cycles the plastic deformation in plane stress conditions vanishes; the macroscopic behaviour of the metal is elastic except around flaws, metallurgical defects or inclusions.

In very high cycle fatigue, the plane stress conditions are not enough for a surface plastic deformation according the Von Mises criteria. The initiation may be located in an internal zone. When the crack initiation site is in the interior, this leads to the formation of one fish eye on the fracture surface, typical of gigacycle fatigue. In this case, the cyclic plastic deformation is related to the stress concentration around a defect: inclusion, porosity, supergrain. Because the probability of occurrence of a flaw is greater within an internal volume than at a surface, the typical initiation in gigacycle fatigue will be most often in the bulk of the metal. During VHCF, the cracks nucleate at microscopic flaws and defects, pore, inclusion, and the crack growth leads to the formation of a structure so-called fish-eye, which could be described as a circular pattern that develops on the fracture surfaces. Ultrasonic fatigue testing 1,4,5 is a powerful tool for evaluating VHCF properties. It uses a high frequency of 20 kHz, which complete 10 9 cycles in a day. The crack propagation in the tested specimen leads to a cyclic plastic dissipation and thus to a source of heat generation. One way of assessing this dissipated energy is to use temperature measurements because part of this dissipation will occur as heat. To this end infrared cameras can be used. [START_REF] Larosa | Thermographic methodology for rapid determination of the fatigue limit of materials and mechanical components[END_REF][START_REF] Fargione | Rapid determination of the fatigue curve by the thermographic method[END_REF] As far we know few studies have been devoted to the computation of the temperature field associated with the propagation of a fatigue crack during VHCF failure. [START_REF] Ranc | Study of thermal effects associated with crack propagation during very high cycle fatigue tests[END_REF][START_REF] Ranc | Thermal effect of plastic dissipation at the crack tip on the stress intensity factor under cyclic loading[END_REF] Our research group has proposed a first simplified analytical approach to model the thermal dissipation during the fish eye propagation, [START_REF] Ranc | Study of thermal effects associated with crack propagation during very high cycle fatigue tests[END_REF] this modification assume that the thermal dissipation is distributed uniformly around the crack tip. The purpose of this work is to enhance this first approach by introducing a mechanical modelling of the cyclic plastic dissipation, without assuming a priori an uniform distribution. We propose a finite element based computational model of the crack propagation during a VHCF failure, in order to compute the temperature field evolution in the tested specimen. The crack propagation is modelled by a classical Paris-Hertzberg crack growth law, and the plastic dissipation per cycle is determined by 3D-elastic plastic finite element modelling of stationary mode I crack under constant amplitude loading. A fraction of this plastic dissipation is converted heat and used as a mobile thermal source to compute the temperature field evolution. The paper is organized as follows: In the first section, the crack growth law and the thermomechanical model are developed. The next section is devoted to the finite element modelling. Finally, the last section the temperature field obtained by the proposed approach is compared with experimental data and the evolution of the plastic dissipation as well as the temperature field around the crack are studied.

M O D E L D E S C R I P T I O N

We consider a cylinder C which has a circular crosssection of radius R c and a height 2L. A small circular crack perpendicular to the cylinder axis lies in the centre of the cylinder, as shown in Fig. 1. The radius of the crack is denoted by a(t) and its eccentricity from the centre of C is denoted by e. C is submitted to a cyclic loading: During each cycle the stress is varied linearly from an initial minimum value σmin to a maximum value σmax and back to the initial value σmin . An important parameter used to characterize the cyclic load is the so-called load ratio, defined as R = σmin / σmax . Because of symmetry conditions only one quarter of the cylinder is modelled (Fig. 2).

Crack growth model

Paris and co-workers [START_REF] Bathias | Gigacycle Fatigue in Mechanical Practice[END_REF] have developed estimation for crack growth life for internal initiation, termed fish-eye using the Paris-Hertzberg crack growth rate law. [START_REF] Hertzberg | Deformation and Fracture Mechanics of Engineering Materials[END_REF] In the present work we use this law to describe the growth of the fatigue crack

da dN = b K eff E √ b 3 , ( 1
)
where b is the Burger's vector modulus, E the elastic modulus and K eff the effective stress intensity factor range.

In this first approach we consider the following simplifications in the history of the crack propagation: the crack is supposed to remain circular during the crack growth process, the crack closure effect is neglected and we use the analytical formula for a circular crack in an infinite media. 11

K eff = K = 2 π σ √ πa. ( 2
)
Let us denotes a 0 the crack radius at t = 0, and by the loading frequency. Hertzberg 10 notes that the threshold corner is such that da dN = b when a = a 0 .

(3)

These assumptions lead to the following expression for the crack radius

a(t) = a 0 1 - b t 2a 0 2 . ( 4
)
The crack grows until reaching the surface of the cylin-

der at t = T f when a(T f ) = R c (1 -e). T f = 2a 0 b 1 - a 0 R c (1 -e) . ( 5
)

Thermal model

We consider a problem of a fatigue crack in an elasticplastic material. The temperature field in the cylinder is evaluated by solving the heat transfer equation

ρC ∂ T ∂t = λ T + Ẇ el + Ẇirr -Ẇs ( 6
)
where ρ is the mass density, C the specific heat, λ the heat conduction coefficient. Ẇ el , Ẇirr , Ẇs are respectively, the thermoelastic structural heating, the plastic heating and the energy storage in the material through the creation of micro-stress fields linked to the development of dislocations. Under fatigue loading the thermoelastic source has the same frequency spectrum as the stress signal and the variation in thermoelatic energy W el vanishes at the end of each loading cycle. [START_REF] Morabito | Analysis of heat sources accompanying the fatigue of 2024 t3 aluminium alloys[END_REF] As we are interested in the temperature variation over a great number of cycles, the temperature may be considered averaged over one complete cycle and the mean temperature rise in the material will be due to the inelastic effects only. Thus the thermoelastic heating Ẇ el can be taken out of Eq. ( 6). The energy storage term Ẇs is usually computed by considering that it represents between 5% and 15% of the plastic heating Ẇirr . [START_REF] Taylor | The latent heat remaining in a metal after cold work[END_REF][START_REF] Simo | Associative coupled thermoplasticity at finite strains: formulation, numerical analysis and implementation[END_REF] The two can be merged in a unique expression written as β Ẇirr , where β is a multiplicative factor, called the Taylor-Quinney factor, which classically is taken between 0.85 and 0.95. Mean plastic energy dissipation per cycle is computed [START_REF] Klingbeil | A total dissipated energy theory of fatigue crack growth in ductile solids[END_REF] where is the frequency of the loading, to obtain a simplified model of the heat transfer equation

Ẇirr

ρC ∂ T ∂t = λ T + β cycle σ : dε p . ( 8
)
It was shown that the characteristic time of the heat convection is higher than the crack propagation time and that the losses due to the heat convection could be neglected during the crack propagation. [START_REF] Ranc | Study of thermal effects associated with crack propagation during very high cycle fatigue tests[END_REF] These considerations lead to the introduction of adiabatic conditions on the surface of the specimen. Moreover, symmetry conditions on the crack plane allow adiabatic condition to be chosen on that plane.

Computation of the plastic energy dissipation

The plastic energy dissipation per cycle is obtained by a 3D elastic-plastic finite element analysis. Three plastic regions exist during the cyclic loading, the forward plastic region which forms during the loading, the reverse plastic region which develops when the load is reversed and the plastic wake which develops when the crack extend in the plastic region (Fig. 3). It was shown that in VHCF regime, [START_REF] Marines-Garcia | Fatigue crack growth from small to large cracks on very high cycle fatigue with fish-eye failures[END_REF] the crack growth is not a significant portion of life in gigacycle fatigue with fish-eye failure (less than 1%). However, the number of cycles involved is still important (≈ 10 5 cycles) and a direct nonlinear computation will lead to prohibitive computational time. In this respect, we propose to follow the approach developed by Klingbeil [START_REF] Klingbeil | A total dissipated energy theory of fatigue crack growth in ductile solids[END_REF] and to compute the energy dissipation per cycle dW irr /dN during a single load cycle on a stationary crack for different radius a i (a in = a 0 < a 1 < • • • < a n = a f i ). To simulate the plastic fatigue regime, and not a monotonic loading, two cycles are computed for each radius a i and only the last cycle is used to evaluate the plastic dissipation. Such an approach neglects the contribution of the actual crack extension during any given load cycle. However, for Paris-Regime crack growth, both the plastic work and the surface energy contributions associated with the actual crack extension in any given cycle are negligible com-pared to the total plastic dissipation. In VHCF regime the imposed load is small compared to the yield strength and the contribution of the plastic wake will be neglected in our work. Finally an elastic-perfectly plastic constitutive model is considered in this work. The yield function is defined by

f (σ ) = σ e -σ Y , ( 9
)
where σ Y is the yield stress and σ e is the Von-Mises equivalent stress defined from the deviatoric stress σ D by

σ e = 3 2 σ D : σ D . ( 10
)
During the cyclic loading, part of the crack surfaces, may be under contact and contact between surface is modelled by a unilateral contact law without friction.

F I N I T E E L E M E N T C O M P U T A T I O N

The total plastic dissipation per cycle is obtained by a 3-D elastic-plastic finite element analysis using Cast3M 17 a general FE computer code for the analysis of structures by the finite element method. This code was developed by the Mechanical Department and Technology (DMT) of the French Police station with Atomic Energy (ECA). The finite element computation is done in two steps. In a first step a set of mechanical computations are done for crack radius a i (a in = a 0 < a 1 < • • • < a n = a f i ). The energy dissipation in the structure during one cycle for each selected a i is denoted (dW irr /dN)(a i ). Each of these computations consists in the analysis of two loading cycles. To obtain an accurate value of the energy dissipation at the crack tip we use linear 3D elements, with a sufficient number of elements in the reversed plastic region. [START_REF] Ellyin | 3d modelling of plasticity induced fatigue crack closure-effect of materials constitutive relations[END_REF] Because of symmetry conditions only one quarter of the cylinder is meshed as shown in Fig. 4. Figure 5 shows the mesh in the plane of the crack for a fixed crack length and Fig. 6 is a zoom of the 3D mesh around the crack tip. The number of time step is set to p for each cycle, a trapezoidal integration rule is used to estimate the value of the dissipation during one cycle on each Gauss point in the mesh.

dW irr dN (a i ) = cycle σ : dε p = p j =1 (σ j +1 + σ j ) : ε p j +1 -ε p j . ( 11
)
In a second step, the thermal problem is solved with an implicit integration scheme. The time discretization (t 0 < t 1 < • • • < t n ) is chosen such that a(t i ) = a i thanks to Eq. ( 4). Where t 0 = 0 designates the initiation of the fish-eye (i.e. a(t 0 ) = a 0 ) and t n = T f designates the instant when the crack reach the specimen surface as defined in Eq. ( 5). The loading is given by the energy dissipation per cycle computed in Eq. (11). The meshes used are the meshes used during the mechanical computation to avoid errors due the data transfer on (dW irr /dN)(a i ) between different meshes.

N U M E R I C A L R E S U L T S

Comparison with measured data

In this section the hypothesis chosen to build the numerical model are evaluated by a comparison with an experimental result. The tested material is a high-strength steel SAE 5120. Fatigue test are performed at ultrasonic fatigue frequency = 20 kHz with a stress ratio of R = 0.1 and stress amplitude σ = 400 MPa, using compressed air of 20 • C to cool the specimen. Both the ultrasonic fatigue machine and the specimen were described in previous pa-per. [START_REF] Wu | An automatic ultrasonic fatigue testing system for studying low crack growth at room and high temperatures[END_REF][START_REF] Bathias | Piezoelectric fatigue testing machines and devices[END_REF] During the loading, an infrared camera is used to measure and record the temperature at the surface of the specimen. The acquisition frequency was 100 Hz with each image consisting of 64 × 80 pixels and a spatial resolution of 0.12 mm/pixel. Two tests have been conducted. In Table 1 are given the number of cycles to failure N f , the crack parameters measured after the failure a inc and e, as well as the radius of the crack at t = 0 which is defined by a 0 = a inc /0.94 as proposed by Paris. [START_REF] Paris | The relationship of effective stress intensity, elastic modulus and burger's vector on fatigue crack growth as associated with fish-eye gigacycle phenomena[END_REF] The final time step is computed from Eq. ( 5) where the Burger vector's modulus is taken such that b = 1.8 × 10 -10 m and the radius of the specimen is supposed to be R c = 1.5 mm. The computations lead to T f = 4.84 s for test 1 and T f = 6.20 s for test 2. Figure 7 shows the evolution of the temperature measured on point P 1 (0, R 1 , 0) (see Fig. 1) dur-ing the formation of the fish-eye (i.e. on the time interval [0, T f ]). The measured temperature at the beginning of the initiation of the fish eye was T(t = 0) = 21.9 • C for test 1 and T(t = 0) = 24 • C for test 2. It can be seen that the plastic work leads to a fast increase of the temperature. The material is approximated as elastic-perfectly plastic with thermomechanical properties E = 200 GPa (elastic modulus), ν = 0.3 (Poisson ratio), σ Y = 1067 MPa (yield stress), ρ = 7800 kg m -3 (mass density), C = 460 J K -1 kg -1 (specific heat), λ = 52 W K -1 m -2 (heat conduction coefficient). The radius of the specimen is supposed to be R c = 1.5 mm and its height is defined by L = 6 mm. The coupled thermomechanical problem defined in the preceding is used to compute the evolution of the temperature on the time interval [0, T f ], where t = 0 denotes the crack initiation and t = T f denotes the failure of the specimen. The initial temperature is assumed to be constant in the specimen and set to the measured data: T(t = 0) = 21.9 • C for test 1 and T(t = 0) = 24 • C for test 2. The comparison between the measured and the computed temperature are plotted in Fig. 8 for test 1, and in Fig. 9 for test 2. These comparisons show a good agreement between the predicted and the measured values. In particular, we observe that, accordingly to the experimental results, the final temperature computed in test 1 is greater than the final temperature computed in test 2. Numerically this fact is linked to the eccentricity of the crack. In test 1 the eccentricity which is smaller than in test 2 leads to a greater final radius of the crack, and thus to a greater plastic work and a greater temperature increase.

Numerical study of the plastic dissipation and of the temperature rise

The following results illustrate the evolution of the plastic dissipation per cycle, and of the temperature near the crack tip. The calculations are done with the geometric and materials values used for test 1. The evolution of the plastic dissipation is computed as a function of the crack length. To study the evolution of the plastic dissipation around the crack, we have divided the crack tip line in regular angular sections (0 = θ 0 < θ 1 < • • • < θ N = 180), and the structure in N elementary volumes V([θ i-1 θ i ]) as shown in Fig. 10. The plastic dissipation per cycle is computed in each volume i-1 θ ]) as a function the length.

dW irr dN ([θ i-1 θ i ]) = V([θi-1 θi ] cycle σ : dε p dV . ( 12
)
The evolution of the plastic dissipation for angular sections of 15 • is shown in Fig. 11. It is important to note that the plastic dissipation is not distributed regularly around the crack tip line as the crack progresses and is mainly concentrated in the two elementary volume which are the nearest to the surface. The contours of the temperature field near the crack tip are shown in Figs 12-14 at different times during the crack propagation. It can be seen in Fig. 12 that the temperature increase is distributed regularly around the crack tip, but is very small. The temperature in the specimen begins to increase when the crack front approaches the edge of the specimen and the stress state is not anymore the stress state of an internal crack in infinite medium Figs 13 and 14).

C O N C L U S I O N

This study has proposed a thermomechanical modelling of the propagation of a fatigue crack in VHCF fatigue regime. The plastic dissipation per cycle has been computed from 3D finite element models of stationary crack in a cylinder under constant amplitude, mode I loading. The temperature rise during fatigue crack propagation is deduced from this plastic dissipation by the resolution of a transient heat conduction problem with a moving source. The numerical results are in a good agreement with experimental test, and show that the plastic dissipation, and hence the thermal source, is concentrated near the surface of the specimen.

N

  O M E N C L A T U R E a inc = average size of the inclusion a 0 = crack radius at t = 0 a(t) = radius of the crack b = Burger's vector modulus C = specific heat e = eccentricity of the crack E = elastic modulus N f = number of cycles to failure R = load ratio R c = radius of the specimen Ẇ el = thermoelastic structural Ẇirr = plastic heating Ẇs = energy storage in the material β = Taylor-Quinney factor K eff = effective stress intensity factor range λ = heat conduction coefficient ν = Poisson ratio ρ = mass density σmin = minimum value of the cyclic stress σmax = maximum value of the cyclic stress σ Y = yield stress σ e = Von-Mises equivalent stress = loading frequency I N T R O D U C T I O N

Fig. 1

 1 Fig. 1 Model description.

Fig. 2

 2 Fig. 2 Model description: one quarter of the cylinder.

Fig. 3

 3 Fig. 3 Plastic regions during cyclic loading.

Fig. 4

 4 Fig. 4 Finite element mesh.

Fig. 5

 5 Fig. 5 Projection of the mesh in the crack plane.

Fig. 6

 6 Fig. 6 Finite element mesh around the crack front.

Fig. 7

 7 Fig.7 Evolution of the temperature during the formation of the fish-eye for an SAE5120. Stress amplitude (400 MPa), load ratio (0.1) and testing frequency (20 kHz).

Fig. 8

 8 Fig. 8 Test 1: Comparison between the modeled and the experimental temperature during the formation of the fish-eye (initial inclusion 10.03 μm, eccentricity 0.79).

Fig. 9

 9 Fig.9 Test 2: Comparison between the modeled and the experimental temperature during the formation of the fish-eye (inclusion 13.83 μm, eccentricity 0.83).

Fig. 10

 10 Fig. 10 Division of the crack tip line in the crack plane.

Fig. 11

 11 Fig.11Evolution of the plastic dissipation as a function of the fish-eye radius.

Fig. 12

 12 Fig. 12 Temperature field in • C at t = 4.5 s.

Table 1

 1 Experimental parameters

	Test	N f (cycles)	a inc (μm)	e	a 0 (μm)	T f (s)
	1	5 . 2 5× 10 7	10.03	0.79	10.67	4.84
	2	4 . 2 4× 10 7	13.80	0.83	14.68	6.20

Acknowledgement

The authors would like to express their gratitude to Dr. Zhi Yong Huang (Sichuan University, China), for the experimental results used to assess the numerical model.


Arnold Sanjay Ramassamy
email: sanjay.ramassamy@ens-lyon.fr

MODULAR PERIODICITY OF THE EULER NUMBERS AND A SEQUENCE BY

For any positive integer q, the sequence of the Euler up/down numbers reduced modulo q was proved to be ultimately periodic by Knuth and Buckholtz. Based on computer simulations, we state for each value of q precise conjectures for the minimal period and for the position at which the sequence starts being periodic. When q is a power of 2, a sequence defined by Arnold appears, and we formulate a conjecture for a simple computation of this sequence.

Introduction

The sequence of Euler up/down numbers (E n ) n≥0 is the sequence with exponential generating series (1)

∞ n=0 E n n!
x n = sec x + tan x.

It is referenced as sequence A000111 in [START_REF] Sloane | The online encyclopedia of integer sequences[END_REF] and its first terms are 1, 1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521, 353792, 2702765, . . .

The numbers E n were shown by André [START_REF] André | Développements de sec x et de tang x[END_REF] to count up/down permutations on n elements (see Section 3). Knuth and Buckholtz [START_REF] Donald | Computation of tangent, Euler, and Bernoulli numbers[END_REF] proved that for any integer q ≥ 1, the sequence (E n mod q) n≥0 is ultimately periodic. For any q ≥ 1 we define :

• s(q) to be the minimum number of terms one needs to delete from the sequence (E n mod q) n≥0 to make it periodic ; • d(q) to be the smallest period of the sequence (E n mod q) n≥s(q) . For example, the sequence (E n mod 3) starts with 1, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, . . . so one might expect to have s(3) = 1 and d(3) = 4. Clearly s(1) = 0 and d(1) = 1. In the remainder of this paper, we formulate precise conjectures for the values of s(q) and d(q) for any q ≥ 2.

Organisation of the paper. In Section 2 we reduce the problem to the case when q is a prime power and we conjecture the values of s(q) and d(q) when q is an odd prime power. In Section 3 we conjecture the values of s(q) and d(q) when q is a power of 2, after having introduced the Entringer numbers and a sequence defined by Arnold describing the 2-adic valuation of the Entringer numbers. In Section 4, we provide a simple construction which conjecturally yields the Arnold sequence.

2. Case when q is not a power of 2

The following lemma implies that it suffices to know the values of s(q) and d(q) when q is a prime power in order to know the values of s(q) and d(q) for any q ≥ 2.

Lemma 1. Fix q ≥ 2 and write its prime number decomposition as

(2) q = k i=1 p α i i ,
where k ≥ 1, p 1 , . . . , p k are distinct prime numbers and α 1 , . . . , α k are positive integers. Then

s(q) = max 1≤i≤k s(p α i i ) (3) d(q) = lcm(d(p α 1 1 ), . . . , d(p α k k )). ( 4
)
The proof is elementary and uses the Chinese remainder theorem. When q is an odd prime power, Knuth and Buckholtz [KB67] found the following : Theorem 2 ([KB67]). Let p be an odd prime number.

(1) If p ≡ 1 mod 4, then

d(p) = p -1.
(2) If p ≡ 3 mod 4, then

d(p) = 2p -2. (3) For any k ≥ 1, s(p k ) ≤ k. (4) For any k ≥ 2, d(p k )|p k-1 d(p).
We conjecture the following for the exact values of s(q) and d(q) when q is an odd prime power : Conjecture 1. Let p be an odd prime number.

(1) For any k ≥ 1,

s(p k ) = k.
(2) For any k ≥ 2,

d(p k ) = p k-1 d(p).
Conjecture 1 is supported by Mathematica simulations done for all odd prime powers q < 1000.

3. Entringer numbers and case when q is a power of 2 Formulating a conjecture analogous to Conjecture 1 for powers of 2 requires to define, following Arnold [START_REF] Vladimir | Bernoulli-Euler updown numbers associated with function singularities, their combinatorics and arithmetics[END_REF], a sequence describing the behavior of the 2-adic valuation of the Entringer numbers.

3.1. The Seidel-Entringer-Arnold triangle. The Entringer numbers are a refined version of the Euler numbers, enumerating some subsets of up/down permutations. For any n ≥ 0, a permutation σ ∈ S n is called up/down if for any 2 ≤ i ≤ n, we have σ(i -1) < σ(i) (resp. σ(i -1) > σ(i)) if i is even (resp. i is odd). André [START_REF] André | Développements de sec x et de tang x[END_REF] showed that the number of up/down permutations on n elements is E n . For any 1 ≤ i ≤ n, the Entringer number e n,i is defined to be the number of up/down permutations σ ∈ S n such that σ(n) = i. The Entringer numbers are usually displayed in a triangular array called the Seidel-Entringer-Arnold triangle, where the numbers (e n,i ) 1≤i≤n appear from left to right on the n-th line (see Figure 1). The Entringer numbers can be computed using the following recurrence formula (see for example [START_REF] Richard | Enumerative combinatorics[END_REF]). For any n ≥ 2 and for any 1 ≤ i ≤ n, we have (5) e n,i = j<i e n-1,j if n is even j≥i e n-1,j if n is odd .

3.2. Arnold's sequence. Replacing each entry of the Seidel-Entringer-Arnold triangle by its 2-adic valuation, we obtain an infinite triangle denoted by T (see Figure 2). 0

∞ 0 0 0 ∞ ∞ 0 1 1 0 0 2 1 ∞ Figure 2.
First five lines of the triangle T of 2-adic valuations of the Entringer numbers.

We read this triangle T diagonal by diagonal, with diagonals parallel to the left boundary. For any i ≥ 1, denote by D i the i-th diagonal of the triangle T parallel to the left boundary. For example D 1 starts with 0, ∞, 0, ∞, 0, . . .. For any i ≥ 1, denote by m i the minimum entry of diagonal D i . Arnold [START_REF] Vladimir | Bernoulli-Euler updown numbers associated with function singularities, their combinatorics and arithmetics[END_REF] observed that the further away one moves from the left boundary, the higher the 2-adic valuation of the Entringer numbers becomes. In particular, he observed (without proof) that the sequence (m i ) i≥1 was weakly increasing to infinity. He defined the following sequence : for any k ≥ 1,

u k := max {i ≥ 1|m i < k} .
In other words, u k is the number of diagonals containing at least one entry that is not zero modulo 2 k . The sequence (u k ) k≥1 is referenced as the sequence A108039 in OEIS [START_REF] Sloane | The online encyclopedia of integer sequences[END_REF] and its first few terms are given in Table 1. Note that the first few terms given by Arnold were incorrect, because the entry 4 appeared four times, whereas it should be appearing only three times. We also remark that we cannot define any sequence analogous to (u k ) when studying the p-adic valuations of the Entringer numbers for odd primes p. Indeed, the p-adic valuation 0 seems to appear in diagonals of arbitrarily high index.

k
3.3. Case when q is a power of 2. Using the sequence (u k ) k≥1 , we formulate the following conjecture for s(q) and d(q) when q is a power of 2 : Conjecture 2. For any k ≥ 1, we have

(6) s(2 k ) = u k .
Furthermore, if k ≥ 1 and k = 2, we have

(7) d(2 k ) = 2 k .
Finally, we have d(4) = 2.

Numerical simulations performed on Mathematica for k ≤ 12 support Conjecture 2.

Construction of Arnold's sequence

In this section we provide a construction which conjecturally yields Arnold's sequence (u k ) k≥1 .

We denote by Z + the set of nonnegative integers and we denote by

Figure 1 .

 1 Figure 1. First five lines of the Seidel-Entringer-Arnold triangle.

+

  the set of all finite sequences of nonnegative integers. We define a map f : S → S, which maps each Z d + to Z 2d + , as follows. Fix x = (x 1 , . . . , x d ) ∈ S. If all the x i 's are equal to x d , we setf (x) = (x d , . . . , x d , 2x d , . . . , 2x d ),where x d and 2x d both appear d times on the right-hand side. Otherwise, defines := max {1 ≤ i ≤ d -1|x i = x d } and set f (x) = (x 1 , . . . , x d , x 1 + x d , . . . , x s-1 + x d , 2x d , . . . , 2x d ),where 2x d appears d -s + 1 times on the right-hand side. For example, we have (8) f ((2, 4, 4, 4)) = (2, 4, 4, 4, 8, 8, 8, 8)

Table 1 .

 1 The first few values of u k .

	1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
	u k 2 4 4 4 8 8 8 8 10 12 12 16 16 16 16 16 18 20

Acknowledgments

The author acknowledges the support of the Fondation Simone et Cino Del Duca.

and (9)

f (2, 4, 4, 4, 8, 8, 8, 8) = (2,4,4,4,8,8,8,8,10,12,12,16,16,16,16,16).

By iterating this function f indefinitely, one produces an infinite sequence :

There exists a unique (infinite) sequence (X k ) k≥1 such that for any k ≥ 1 and for any n ≥ log 2 (k/d), X k is the k-th term of the finite sequence f n (x). This infinite sequence is called the f -transform of x. The lemma follows from the observation that for any ≥ 1 and for any y ∈ Z + , y and f (y) have the same first terms.

We can now formulate a conjecture about the construction of the sequence (u k ) k≥1 : Conjecture 3. Arnold's sequence (u k ) k≥1 is the f -transform of the quadruple (2, 4, 4, 4). Conjecture 3 is supported by the estimation on Mathematica of u k for every k ≤ 512.


David Ryckelynck
email: david.ryckelynck@mines-paristech.fr

Laurent Gallimard

Samuel Jules

Estimation of the validity domain of hyper-reduction approximations in generalized standard elastoviscoplasticity

Keywords: Reduced integration domain, Hyper-reduction, POD, Incremental variational principle

Background:

We propose an a posteriori estimator of the error of hyper-reduced predictions for elastoviscoplastic problems. For a given fixed mesh, this error estimator aims to forecast the validity domain in the parameter space, of hyper-reduction approximations. This error estimator evaluates if the simulation outputs generated by the hyper-reduced model represent a convenient approximation of the outputs that the finite element simulation would have predicted. We do not account for the approximation error related to the finite element approximation upon which the hyper-reduction approximation is introduced. Methods: We restrict our attention to generalized standard materials. Upon use of incremental variational principles, we propose an error in constitutive relation. This error is split into three terms including a tailored norm of the hyper-reduction approximation error. This error norm is defined by using the convexity of an incremental potential introduced to state the constitutive equations. The second term of the a posteriori error is related to the stress recovery technique that generates stresses fulfilling the finite element equilibrium equations. The last term is a coupling term between the hyperreduction approximation error at each time step and the errors committed before this time step. Unfortunately, this last term prevents error certification. In this paper, we restrict our attention to outputs extracted by a Lipschitz function of the displacements. Results: In the proposed numerical examples, we show very good preliminary results in predicting the validity domain of hyper-reduction approximations. The average computational time of the predictions obtained by hyper reduction, is accelerated by a factor of 6 compared to that of finite element simulations. This speed-up incorporates the computational time devoted to the error estimation. Conclusions: The numerical implementation of the proposed error estimator is straightforward. It does not require the computation of the incremental potential. In the numerical results, the estimated validity domain of hyper-reduced approximations is inside the reference validity domain. This paper is a first attempt for a posteriori error estimation of hyper-reduction approximations.

Background

Reduced-Order models aim at reducing the computational time required to obtain solutions of Partial Differential Equations (PDE) which are physically-based and parameter dependent. They reduce the computational complexity of optimization procedures or parametric analyses [START_REF] Ganapathysubramanian | A non-linear dimension reduction methodology for generating data-driven stochastic input models[END_REF][START_REF] Balima | Model reduction for heat conduction with radiative boundary conditions using the modal identification method[END_REF][START_REF] Daescu | Efficiency of a pod-based reduced second-order adjoint model in 4d-var data assimilation[END_REF].

Reduced-Order Models are very useful to model complex mechanical phenomena, when parametric studies are mandatory to setup convenient nonlinear constitutive equations or boundary conditions. In such situations, many open questions about physical assumptions arise and the scientist's intuition alone cannot guarantee that convenient simplified numerical approximations are chosen. Here, we consider hyper-reduced (HR) simulations, involving both a reduced-basis approximation and reduced-coordinate estimation by using a mesh restricted to Reduced Integration Domain (RID) [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF]. The introduction of the RID is crucial for elastoviscoplastic models, because in many practical cases, no speed-up is achieved if the mesh is not restricted to the RID. In such a framework, error estimation provides a valuable algorithmic approach to check if hyperreduced simulations have been performed in the validity domain of the hyper-reduced model, and if they allow a physical understanding of the simulation outputs.

HR simulations may not be performed in faster time than predictions using metamodeling, response surface methods, or virtual charts [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer[END_REF]. But, unlike these methods, Reduced-Order Models involve equations related to physical principles (i.e. the balance of momentum), physical constitutive equations and data, validated by experiments. Such an approach enables certification of outputs of interest [START_REF] Cuaong | Certified real-time solutin of parametrized partial differential equations[END_REF], or error estimation of predictions obtained by using Reduced-Order Models [START_REF] Cuaong | Certified real-time solutin of parametrized partial differential equations[END_REF][START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF][START_REF] Kerfriden | Certification of projection-based reduced order modelling in computational homogenisation by the constitutive relation error[END_REF][START_REF] De Almeida | A basis for bounding the errors of proper generalised decomposition solutions in solid mechanics[END_REF]. Moreover, error estimation is mandatory in many model reduction of parametric Partial Differential Equations, when one needs to evaluate, in the parameter space, a trust region or a validity domain related to Reduced-Basis accuracy [START_REF] Arian | Trust-region proper orthogonal decomposition for flow control[END_REF][START_REF] Michel | Optimal control of the cylinder wake in the laminar regime by trust-region methods and pod reduced-order models[END_REF][START_REF] Bui | Pod-isat: An efficient pod-based surrogate approach with adaptive tabulation and fidelity regions for parametrized steady-state pde discrete solutions[END_REF]. Error estimation is also mandatory to perform a priori model reduction by using greedy algorithms [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF][START_REF] Bris | Results and questions on a nonlinear approximation approach for solving high-dimensional partial differential equations[END_REF] or the A Priori Hyper-reduction method [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF].

Unfortunately, in nonlinear mechanics of materials, parametric PDEs are rarely affine in parameter, and it is quite difficult to reduce equations and data to an approximate affine form as proposed in [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF] by using the Empirical Interpolation Method [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Nguyen | A 'best points' interpolation method for efficient approximation of parametrized functions[END_REF]. The novelty of this paper is the use of an incremental variational principle [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF][START_REF] Miehe | Strain-driven homogenization of inelastic micro-structures and composites based on an incremental variational formulation[END_REF][START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I. incremental variational principles[END_REF] to develop an a posteriori estimator of errors for HR predictions in case of generalized standard elastoviscoplastic models. As usual, the finite element (FE) solution fulfills incremental equations in time, but no differential equations related to the continuous time variable. The incremental variational principle enables to preserve formulations discrete in time when comparing the HR solution to the FE solution. This approach enables to state a relationship between a tailored norm of the hyper-reduction approximation error and the a posteriori estimator of this error. Here, we term hyper-reduction approximation error (HRAE) the difference between the FE solution of the PDE and its estimation by the hyper-reduced model. The reference solution is the FE solution, not the exact solution of the PDE. In the sequel, the proposed tailored norm of the displacement is denoted by u u , where u is the predicted displacement over the domain . The error estimator is denoted by η. We follow the theory of constitutive relation error (CRE) proposed by P. Ladevèze [START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF]. This kind of a posteriori error has been developed for a large set of constitutive equations [START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF][START_REF] Gallimard | Error estimation and adaptivity in elastoplasticity[END_REF][START_REF] Ladevèze | A posteriori constitutive relation error estimators for nonlinear finite element analysis and adaptive control[END_REF][START_REF] Ladevèze | Strict upper bounds of the error in calculated outputs of interest for plasticity problems[END_REF][START_REF] Bouclier | Real-time validation of mechanical models coupling pgd and constitutive relation error[END_REF]. Although constitutive relation error (CRE) can handle all approximation errors due to time discretization, FE approximation, and eventually a separated representation of the displacements [START_REF] Pelle | An efficient adaptive strategy to master the global quality of viscoplastic analysis[END_REF], we restrict our attention to approximation errors due to the additional assumptions incorporated in the numerical model by the hyper-reduced formulation. Here, the CRE η requires the construction of stresses fulfilling FE equilibrium equations.

In the numerical part of this paper, we apply the error estimation to forecast the validity domain of hyper-reduced predictions of simulation outputs. Adjoint sensitivity analyses have been introduced for certification of outputs in [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF][START_REF] Paraschivoiu | A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations[END_REF]. In case of time dependent problems, the adjoint equation or the mirror problem [START_REF] Ladevèze | Strict upper bounds of the error in calculated outputs of interest for plasticity problems[END_REF] is solved backward in time from T to 0. Here, we circumvent such a complex numerical solution by restricting our attention to Lipschitz functions of the displacements for simulation outputs, similarly to [29]. The output function is denoted by s(u) ∈ R N s . The Lipschitz continuity of the output s(u) reads:

∃ c o ∈ R + s(u) -s(v) 2 F c o u -v 2 u , u ∈ V h , v ∈ V h , ( 1 )
where V h is the convenient Hilbert space related to the FE approximation of the PDE, • F is the Frobenius norm on R N s . The elastoviscoplastic problem is both discretized in space and in time by the finite element method and the forward Euler scheme. The solution of this FE problem at discrete time t n ∈ {t o , . . . , t N t } is denoted by u n FE :

u n FE (x; μ) = N ξ i=1 ξ i (x) q n i (μ) x ∈ , μ ∈ D, V h = span(ξ 1 , . . . , ξ N ξ ) (2)
where

(ξ i ) N ξ
i=1 are the FE shape functions, x the space coordinate, μ is a vector of parameters and D is a given parameter space.The hyper-reduction approximation of this FE solution is denoted by u n HR :

u n HR (x; μ) = N ψ k=1 ψ k (x) γ n k (μ) x ∈ , μ ∈ D, ψ k ∈ V h (3)
where

(ψ k ) N ψ k=1 is a given reduced-basis, (γ n k ) N ψ
k=1 are the reduced coordinates evaluated by using equilibrium equations set-up on the RID denoted by Z , Z ⊂ , by following the hyper-reduction method [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF]. The HRAE is denoted by e n , e n = u n FEu n HR . Because of the time discretization, we cannot consider the approximation error e n separately for each time index n. We must account for the error at all discrete times simultaneously. Therefore, the upper bound of interest related to the Lipschitz continuity reads:

∃ c o ∈ R + N t n=1 s(u n HR + e n ) -s(u n HR ) 2 F c o N t n=1 e n 2 u , e n ∈ V h , u n HR ∈ V h (4)
The reference validity domain of the hyper-reduced model is denoted by D V . Its definition reads:

D V = μ ∈ D | N t n=1 s u n HR + e n -s(u n HR ) 2 F D ( 5
)
where D is a given tolerance. The proposed numerical error estimation aims to forecast a reliable approximation of D V by the recourse to hyper-reduced simulations, a pos-teriori estimations of errors and a FE simulation used to generate reduced-bases. The approximate validity domain is denoted by D V :

D V = μ ∈ D | c η η (u n HR ) N t n=1 , ( σ n ) N t n=1 D (6)
where c η is a substitute for the Lipschitz constant c o , η is a substitute for the tailored norm of the HRAE, (u n HR ) N t n=1 is the sequence of hyper-reduced predictions over the full time interval and ( σ n ) N t n=1 is a sequence of stresses fulfilling the FE equilibrium equation at each discrete time t n , where σ n is an affine function of the stress estimated from u n HR . In the sequel, both sequences (u n HR ) N t n=1 and ( σ n ) N t n=1 are denoted by u HR and σ respectively. Various Reduced-Basis construction methods can be incorporated into the hyperreduced model. For the sake of simplicity, in the following numerical examples, the Reduced-Basis is provided by a Proper Orthogonal Decomposition [START_REF] Karhunen | Über lineare methoden in der wahrscheinlichkeitsrechnung[END_REF][START_REF] Loève | Probability Theory[END_REF][32] of the FE solution related to a given parameter value μ 1 ∈ D. The POD modes are orthonormal solutions of the following mimization problem:

(ψ k ) N ψ k=1 = arg min (ψ k ) N ψ k=1 ∈V h N t n=1 u n FE (•; μ 1 ) -ψ k ψ k , u n FE (•; μ 1 ) 2 L 2 ( ) (7)
where •, • is the L 2 scalar product defined on and • L 2 ( ) is the related norm. The hyper-reduced model is set-up in order to provide an accurate approximate solution for μ = μ 1 :

c η η u HR •; μ 1 , σ •; μ 1 HR , HR D ( 8
)
where HR is a given tolerance. Hence μ 1 ∈ D V . Moreover c η fulfills the following constraint:

N t n=1 s u n FE •; μ 1 -s u n HR •; μ 1 2 F c η η u HR •; μ 1 , σ •; μ 1 (9)
In this paper μ 1 is arbitrary chosen. Nevertheless, the knowledge of an approximate validity domain aims to validate this choice or to adapt it. In the sequel, for clarity, we do not mention the vector of parameters μ in the notations.

Methods

This section is organized as follows. The generalized standard constitutive equations are presented first. They are based on an incremental variational principle proposed in [START_REF] Radovitzky | Error estimation and adaptive meshing in strongly nonlinear dynamic problems[END_REF]. Then, we introduce the hyper-reduced incremental problem to be solved. The following section is devoted to the construction of the FE-equilibrated stress within the framework of the hyper-reduced modeling. The next sections introduce the error estimator, the tailored norm related to the CRE and the partition of the CRE. We finish by remarks on the numerical implementation of the proposed approach.

Incremental potential for the constitutive equations

The constitutive laws are described by using an incremental potential in the framework of the irreversible thermodynamic processes. A priori error estimators and incremental variational formulations were introduced in [START_REF] Radovitzky | Error estimation and adaptive meshing in strongly nonlinear dynamic problems[END_REF] for mechanical problems of bodies undergoing large dynamic deformations. Extensions of this approach were proposed in [START_REF] Miehe | Strain-driven homogenization of inelastic micro-structures and composites based on an incremental variational formulation[END_REF][START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I. incremental variational principles[END_REF] for effective response predictions of heterogeneous materials. The strain history is taken into account by using internal variables denoted by α. These variables are the lump sum of the history of material changes. This approach has its roots in the works by Biot [START_REF] Biot | Mechanics of Incremental Deformations[END_REF], Ziegler [START_REF] Ziegler | Some Extremum Principles in Irreversible Thermodynamics with Applications to Continuum Mechanics[END_REF], Germain [START_REF] Germain | Continuum thermodynamics[END_REF] or Halpen and Nguyen [START_REF] Halphen | Generalized standard materials[END_REF] and has proven its ability to cover a broad spectrum of models in viscoelasticity, viscoplasticity, plasticity and also continuum damage mechanics. The FE solution is approximated by an hyperreduced predictions denoted by u n HR , α n HR , σ n HR respectively for the displacements, the internal variables and the Cauchy stress, at discrete time t n . For the sake of simplicity, we denote by ε n HR the strain tensor ε(u n HR ), in the framework of the infinitesimal strain theory.

According to the framework of the incremental variational formulations, the constitutive law is defined by a condensed incremental potential, denoted by w (ε n HR ), such that:

σ n HR = ∂w ∂ε (ε n HR ), n = 1, . . . , N t (10)
In practice, we do not implement the computation of the partial derivative of w with respect to the component of the strain tensor. Nevertheless, Equation ( 10) is fulfilled up to the computer precision. The internal variables are solutions of the following equation [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I. incremental variational principles[END_REF]:

α n HR = arg Inf α J(ε n HR , α ), (11)
J(ε, α) = w(ε, α) + t n -t n-1 ϕ α -α n-1 HR (x) t n -t n-1 , ( 12
)
w (ε) = Inf α J(ε, α ), ( 13
)
where w(ε, α) is the free energy of the material, and ϕ( α) is its dissipation potential [START_REF] Germain | Continuum thermodynamics[END_REF].

The two potentials w and ϕ are convex functions of their arguments (ε, α) and α respectively, according to the theory of generalized standard materials [START_REF] Germain | Continuum thermodynamics[END_REF][START_REF] Halphen | Generalized standard materials[END_REF]. Examples of constitutive laws can be found in [START_REF] Lemaitre | Mecanique des Materiaux Solides. 1st edn. Dunod, English version[END_REF]. A detailed example is given in the last section of this paper.

The convexity of w is proved in [START_REF] Lahellec | On the effective behavior of nonlinear inelastic composites: I. incremental variational principles[END_REF] under the assumption that w and ϕ are convex functions. In the sequel, the explicit knowledge of the condensed incremental potential is not required for the mathematical formulation of the error estimator.

By recourse to the Legendre transformation, the dual potential of w , denoted by w , reads:

w ( σ ) = Sup ε ε : σ -w ε ( 14
)
We restrict our attention to convex functions w , hence:

w ( ε) = Sup σ ε : σ -w (σ ) (15)
Therefore:

w ( ε) + w ( σ ) -ε : σ ≥ 0 ∀ ε, σ (16)
The definition of the partial derivatives ∂w ∂ε and ∂w ∂σ is extended to fulfill the following equations:

σ ε = arg Sup σ ε : σ -w (σ ) ⇔ ε = ∂w ∂σ ( σ ε ) (17) ε σ = arg Sup ε ε : σ -w (ε ) ⇔ σ = ∂w ∂ε ( ε σ ) (18)
Hence, by construction, σ n HR and ε n HR fulfill the following equation:

w (ε n HR ) + w (σ n HR ) -ε n HR : σ n HR = 0 (19)

Hyper-reduced setting

The continuous medium is occupying a domain . The boundary ∂ of is denoted by ∂ U ∪ ∂ F . On ∂ U , there is the Dirichlet condition u = u c . On ∂ F , there is a given force field F n . The FE displacement field belongs to a function space u c + V h . The reduced subspace is denoted by V ROM :

V ROM = span(ψ 1 , . . . , ψ N ψ ) ⊂ V h ( 20
)
We denote by V the matrix form of empirical modes defined on the FE basis:

ψ k (x) = N ξ i=1 ξ i (x) V ik , k = 1, . . . , N ψ (21)
By following the hyper-reduction method [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF], we generate the RID Z by using the empirical modes (ψ k ) N ψ k=1 and the strain modes (ε(ψ k )) N ψ k=1 . When the RID is known, we determine the set of available FE residual entries when the stress estimation is restricted to Z . This set is denoted by Z such that:

Z = {i ∈ {1, . . . , N ξ }, \ Z ξ 2 i d = 0} (22)
Therefore, we can define truncated test functions insuring a weak form of the balance of momentum restricted to the domain Z , since these test functions are set to zero over \ Z . These truncated test functions are denoted by (ψ Z k )

N ψ k=1 such that:

ψ Z k (x) = i∈Z ξ i (x) V ik , k = 1, . . . , N ψ ( 23
)
The statement of the hyper-reduced incremental problem is the following. Given a parameter vector μ, find an estimation of the reduced coordinates γ n such that u n HR ∈ u c + V ROM fulfills the constitutive equations and the principle of virtual work:

u n HR (x) = u c + N ψ k=1 ψ k (x) γ n k x ∈ (24) Z ε(ψ Z k ) : σ n HRZ d = ∂ F ∩∂ Z ψ Z k . F n d ∀ k ∈ {1, . . . , N ψ } ( 25
)
α n HR = arg Inf α J(ε(u n HR ), α ) ∀x ∈ Z ( 26
)
σ n HRZ = ∂w ∂ε (ε(u n HR )) ∀x ∈ Z ( 27
)
where : is the contracted product for second-order tensors. We must emphasize the fact that Equation [START_REF] Ladevèze | A posteriori constitutive relation error estimators for nonlinear finite element analysis and adaptive control[END_REF] gives access to the displacement in all the domain although the equi-librium is set only on Z . Hence, when the hyper-reduced solution is known, the stress σ n HR can be estimated on the full domain by using Equation [START_REF] Kerfriden | Certification of projection-based reduced order modelling in computational homogenisation by the constitutive relation error[END_REF].

The FE equations are obtained by substituting in equations ( 24) to ( 26), the following items:

• the shape functions (ξ i )

N ξ i=1 for (ψ k ) N ψ k=1 , • the test functions (ξ i ) N ξ i=1 for (ψ Z k ) N ψ k=1 , • the full domain for Z , • the variables ((q i ) N ξ i=1 , u n FE , σ n FE , α n FE ) for (γ k ) N ψ k=1 , u n HR , σ n HR , α n HR .
However, in Equation ( 27), the incremental potential related to the hyper-reduced prediction is not the incremental potential of the FE prediction, because α n-1 FE and α n-1 HR may differ. So, σ n FE fulfills the constitutive equation, but it may not be a derivative of w . Therefore, we introduce a correction term in stress, denoted by δσ n , such that:

σ n FE = ∂w ∂ε (ε n FE ) + δσ n ( 28
)
where ε n FE = ε(u n FE ) and δσ n account for the variation of the internal variables due to the difference between ε i FE n-1 i=1 and ε i HR n-1 i=1 at time instants before t n . The convexity of w and φ insures that the FE solution is unique, if no rigid mode is available.

Dual reduced-subspace

In this section, we introduce the technique of construction of the stress fields σ n fulfilling the finite element equilibrium equation at each discrete time t n . These stress fields are used in the next sections to build a CRE. The FE equilibrium equation is a linear equation upon stresses denoted by σ n . It defines an affine space such that σ n ∈ σ n N + S h , where σ n N is a particular solution of Neumann boundary conditions (i.e. in the linear elastic case for instance), and S h is the following vector space:

S h = σ ∈ L 2 ( ) | ε(u ) : σ d = 0 ∀u ∈ V h ( 29
)
The following linear elastic problem gives us σ n N :

u n N ∈ V h (30) ε(u ) : σ n N d = ∂ F u . F n d ∀ u ∈ V h ( 31
)
σ n N = E o ε(u n N ) ∀x ∈ (32)
where E o ∈ R + is an abritary constant. By the recourse to the stress predicted by the FE simulation for μ = μ 1 , we obtain the snapshots of stress fields

σ n FE (•; μ 1 ) -σ n N N t
n=1 that span a subspace of S h :

σ n FE •; μ 1 -σ n N ∈ S h , n = 1, . . . , N t ( 33
)
Similarly to the approach proposed in [START_REF] Kerfriden | Certification of projection-based reduced order modelling in computational homogenisation by the constitutive relation error[END_REF] for elasticity, we introduce a dual reducedsubspace denoted by S ROM ⊂ S h . S ROM is generated by the usual POD applied to

σ n FE (•; μ 1 ) -σ n N N t n=1
. Its dimension is denoted by N σ ψ , and it is such that N σ ψ ≤ N t .

Hence, the hyper-reduced predictions of the stress tensor can be projected into the space of admissible stresses fulfilling the FE equilibrium equation. This projection reads:

σ n = σ n N + arg min σ ∈S ROM σ n HR -σ n N -σ L 2 ( ) (34)
The above minimization problem is a global L 2 projection of the stress σ n HR onto the reduced basis that span S ROM . The dual basis being a POD basis, it is orthonormal with respect to the L 2 scalar product. So the computational complexity of the stresses projection scales linearly with the number of Gauss points of the mesh and N σ ψ . In practice, this computational complexity is negligible compared to the complexity of the evaluation of σ n HR , by using the constitutive equations.

The constitutive relation error and its partition

The reference solution being incremental in time, no continuous formulation in time is considered for the error estimator. The error estimation proposed in [START_REF] Radovitzky | Error estimation and adaptive meshing in strongly nonlinear dynamic problems[END_REF] for incremental variational formulation, is related to an asymptotic convergence assumption in order to get an upper bound of the approximation error related to the FE discretization. This upper bound depends on a constant related of the weak form of the partial differential equations.

Here, we propose to apply the Constitutive Relation Error method proposed in [START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF] to estimate the constant c η without assuming an asymptotic convergence of the HRAE. Following the method proposed in [START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF], a constitutive relation error, denoted by η, can be introduced as follows:

η(u HR , σ ) = N t n=1 w (ε(u n HR )) + w ( σ n ) -ε(u n HR ) : σ n d ( 35
)
where σ n is defined by Equation [START_REF] Biot | Mechanics of Incremental Deformations[END_REF]. The following properties hold:

η( u, σ ) ≥ 0 ∀ u ∈ u c + V h , ∀ σ (36) η(u HR , σ ) = 0 ⇔ σ n = ∂w ∂ε (ε(u n HR )) ∀x ∈ ∀ n ∈ {1, . . . , N t } (37) η(u HR , σ ) = 0 ⇒ e n = 0 ∀x ∈ ∀n ∈ {1, . . . , N t } (38)
The first property (36) comes from the Legendre transformation [START_REF] Bui | Pod-isat: An efficient pod-based surrogate approach with adaptive tabulation and fidelity regions for parametrized steady-state pde discrete solutions[END_REF]. The property [START_REF] Halphen | Generalized standard materials[END_REF] comes from Equation [START_REF] Ortiz | The variational formulation of viscoplastic constitutive updates[END_REF]. The proof of Property ( 38) is the following. If η(u n HR , σ n ) = 0 ∀ n ∈ {1, . . . , N t } then (u n HR , σ n , α n HR ) fulfills the constitutive equations, the Dirichlet conditions and the FE equilibrium equation. As (u n HR , σ n , α n HR ) fulfills all the equations of the FE problem, it is a solution of the FE problem and e n = 0.

We propose for e n u an Hilbert norm parametrized by the approximate solution u n HR and the exact solution u n FE . This parametrized norm is defined by:

e n 2 u = ε(e n ) : G(ε(u n HR ), ε(u n FE )) : ε(e n ) d ( 39
)
where ε(•) is the symmetric part of the gradient and G a symmetric positive-definite fourth-order tensor. If the identity tensor is substituted for G, we obtain a usual H 1 ( ) norm. We assume that there is no rigid mode, neither in the FE solution nor in the reduced basis.

The following mathematical developments aim to establish a relation between η(u HR , ∂w ∂ε (ε FE )) and the tailored norm e n u . We choose G by the recourse to the convexity of the incremental potential w . According to the Legendre transformation, we have:

w ∂w ∂ε ε n FE = ε n FE : ∂w ∂ε ε n FE -w ε n FE , n = 1, . . . , N t , (40)
Hence, we can remove the contribution of the dual potential w * in η u HR , ∂w ∂ε (ε FE ) :

η u HR , ∂w ∂ε (ε FE ) = N t n=1 w (ε n HR ) -w (ε n FE ) - ∂w ∂ε (ε n FE ) : (ε n HR -ε n FE ) d ( 41
)
Let us define the scalar function f (λ):

λ ∈[ 0, 1] , f (λ) = η β(λ), ∂w ∂ε (ε FE ) , (42) w
.r.t β n (λ) = u n FE + λ u n HR -u n FE , n = 1, . . . , N t (43)
We assume that f is C 2 on [ 0, 1]. The application of the Taylor Lagrange formula at order 2 leads to:

∃ λ c ∈] 0, 1[ s.t. f (1) = f (0) + f (0) + 1 2 f (λ c ) (44)
The function evaluation and its derivatives read:

f (0) = 0 (45) f (1) = η u HR , ∂w ∂ε (ε FE ) (46) f (0) = 0 (47) f (λ) = N t n=1 ε n HR -ε n FE : ∂ 2 w ∂ε 2 ε β n (λ) : ε n HR -ε n FE d (48)
Therefore, by choosing:

G ε n HR , ε n FE = 1 2 ∂ 2 w ∂ε 2 ε n FE + λ c ε n HR -ε n FE (49)
we obtain the following property:

N t n=1 e n 2 u = η u HR , ∂w ∂ε (ε FE ) ( 50
)
This property is an intermediate result before establishing the relationship between e n 2 u and the error estimator. The incremental potential being strictly convex, G is definite positive. A schematic view of N t n=1 e n 2 u is shown on Figure 1. Similarly, a norm on stresses can be defined to measure the distance between the admissible stress σ n and σ n FE -δσ n , by the recourse to the dual potential w such that:

σ n -σ n FE + δσ n 2 σ = σ n -σ n FE + δσ n : H σ n , σ n FE -δσ n : σ n -σ n FE + δσ n d (51)
(ε ) = Nt n=1 w ε n FE + ∂w ∂ε ε n FE : ε n -ε n FE d .
where H is a positive symmetric fourth order tensor. To establish the relation between H and the Hessian matrix of the potential w , we consider η(u FE , σ ). Intermediate stresses (S n (λ)) N t n=1 are introduced such that:

S n (λ) = σ n FE -δσ n + λ ( σ n -σ n FE + δσ n ), n = 1, . . . , N t (52)
By the recourse to the Taylor Lagrange formula, we obtain the following property:

λ c ∈] 0, 1[ , s.t N t n=1 σ n -σ n FE + δσ n 2 σ = η(u FE , σ ) (53) w.r.t H σ n , σ n FE -δσ n = 1 2 ∂ 2 w ∂σ 2 (S n (λ c ))
Property. The proposed constitutive relation error fulfills the following partition:

η(u HR , σ ) = N t n=1 e n 2 u + N t n=1 σ n -σ n FE + δσ n 2 σ + N t n=1 ε(e n ) : δσ n d (54)
The last term of the sum is a coupling term between the error e n and the HRAE committed before the discrete time t n . Since we can't certify that this term is positive, η is not an upper bound of the HRAE.

Thanks to the definition of δσ by Equation ( 28), we have the following properties:

η(u FE , σ FE -δσ ) = 0 (55) η(u HR , σ FE -δσ ) = N t n=1 e n 2 u ( 56
)
We start proving the partition property by considering the following sum of terms:

A = η(u HR , σ ) -η(u HR , σ FE -δσ ) -η(u FE , σ ) + η(u FE , σ FE -δσ )
In A, the contributions of w and w are canceled when considering the definition of η, and we obtain:

A = N t n=1 -ε n HR : σ n + ε n HR : σ n FE -δσ n + ε n FE : σ n -ε n FE : σ n FE -δσ n d
Then:

A = N t n=1 -ε n HR : σ n -σ n FE + δσ n + ε n FE : σ n -σ n FE + δσ n d = N t n=1 ε n FE -ε n HR : σ n -σ n FE d + N t n=1 ε n FE -ε n HR : δσ d , where ε n FE -ε n HR = ε(e n ).
Therefore, according to properties (53), ( 55) and ( 56), the following partition is fulfilled:

η (u HR , σ ) = N t n=1 e n 2 u + N t n=1 σ n -σ n FE + δσ n 2 σ + N t n=1 ε(e n ) : δσ n d + N t n=1 ε(e n ) : σ n -σ n FE d
Since e n ∈ V h and σ n FE -σ n ∈ S h , then:

ε(e n ) : σ n -σ n FE d = 0 ∀n
This ends the proof.

Property. If e n = 0, for n ∈ {1, . . . , N t } then u n HR = u n FE , α n HR = α n FE , δσ n = 0 and:

η(u FE , σ FE ) = N t n=1 σ n FE -σ n FE 2 σ
(57)

σ n FE = σ n N + arg min σ ∈S ROM σ n FE -σ n N -σ L 2 ( ) ( 58
)
In the general case, the closer σ n to σ n FE -δσ n the better the error estimation. The proposed error estimator incorporates errors related to the projection of the stress onto the dual reduced basis.

Numerical implementation of the CRE

The numerical implementation of the incremental potential and its dual is not required for the error estimation. Let's consider the following function of a scalar coordinate λ:

λ ∈[ 0, 1] , g(λ) = w ε n HR + w (τ (λ)) -ε n HR : τ (λ), ( 59
)
τ (λ) = σ n HR + λ σ n -σ n HR ( 60
)
Then:

g(0) = 0 (61) η(u HR , σ ) = N t n=1 g(1) d (62) d g d λ (λ) = ∂w ∂σ (τ (λ)) -ε n HR : σ n -σ n HR (63)
Therefore, the following property holds:

η(u HR , σ ) = N t n=1 1 0 ( ε n (λ) -ε n HR ) : σ n -σ n HR dλ d (64) w.r.t. ε n (λ) = ∂w ∂σ (τ (λ)), τ (λ) = σ n HR + λ σ n -σ n HR (65)
Here, ε n (λ) is the strain tensor given by the constitutive equation upon the stress τ (λ) at time t n and the internal variables α n-1 HR at time t n-1 . Hence, the numerical estimation of ε n (λ) can be performed by the usual implementation of the constitutive equations of generalized standard materials [START_REF] Besson | Non-linear mechanics of materials[END_REF][START_REF] Fritzen | Nonlinear homogenization using the nonuniform transformation field analysis[END_REF][START_REF] Fritzen | Reduced basis hybrid computational homogenization based on a mixed incremental formulation[END_REF].

In the following numerical simulations, we have estimated the integral on λ by the value of first derivative of g for λ = 1.

Numerical estimation of the constant c η

In practice, we have more modes than necessary to achieve accurate predictions for μ = μ 1 , because the POD of (u n FE (•, μ 1 )) N t n=1 is accurate up to the numerical precision. The total number of available modes is denoted by N ψ such that N ψ ≤ min(N ξ , N t ) and N ψ ≤ N ψ . This gives access to several hyper-reduced models by restricting the number of modes involved in the reduced basis related to the displacements. The constant c η is setup by using these approximate hyper-reduced solutions and the known stresses predicted by the FE model:

c η = max N ψ ∈{1,...,N ψ } N t n=1 s u n FE •; μ 1 -s u n HR •; μ 1 2 F min η u HR •; μ 1 , σ FE •; μ 1 , η u HR •; μ 1 , σ •; μ 1 (66)
Hence the constraint (9) is fulfilled. Let us denote by N c ψ the number of displacement modes for which the maximum in Equation (66) is reached. In our opinion, if we expect that the error estimator behaves like an upper bound, we should not take N c ψ modes to generate the final HR model. In the following example, we are setting N ψ = N c ψ + 1.

Results and discussion

We apply the estimation of the validity domain to a beam composed with a layer of tin alloy (SnX) on a layer of copper (Cu). This is a cantilever beam shown in Figure 2. The thickness of the tin is 0.155 mm, the thickness of the copper in 0.111 mm, the width and the length of the beam are respectively 2. mm and 12. mm. These are the dimensions of an experimental specimen used at the Centre des Matériaux. The simulation output is the difference between an experimental transverse displacement at point A, shown in Figure 2, and its prediction by the mechanical model:

s(u n HR ) = y • u n HR (A; μ) -u n exp ( 67
)
where y is the transverse axis and u n exp is a given experimental measurement of the transverse displacement of the point A. The copper is assumed to be elastic and isotropic. Its Young modulus is 63000. MPa and its Poisson coefficient is 0.3. The tin alloy is isotropic and elastoviscoplastic with a linear kinematic hardening. The potentials of the tin alloy read:

α = (ε p , α) (68) w(ε, α) = 1 2 E 1 + ν (ε -ε p ) 2 + ν E (1 -2ν) (1 + ν) Tr(ε -ε p ) 2 + 1 3 C α : α (69) ϕ( α) = 2 3 εp : εp 1 2 R o + h 2 3 ( εp : εp ) 1 2m - h m + 1 2 3 εp : εp m+1 2m + 1 δ((ε p -˙ α) : (ε p -˙ α)) + 1 δ(|Tr(ε p )|) (70) δ(0) = + ∞ (
F n = -F max sin π n N t y, F max = 0.04 N, N t = 5 (73)
The FE simulation used to generate the reduced basis is performed on an arbitrary sampling point μ 1 = (20. MPa, 30000. MPa). The POD of the FE solution gives N ψ = 9 displacement modes with non-zero eigenvalue. It also provides N k=1 the statically admissible modes related to the POD of (σ n FE (•; μ 1 ) -σ n N ) N t n=1 . Some of these modes are shown in Figure 3. Figure 4 shows that the proposed error estimator incorporates errors related to the projection of the stress onto the dual reduced basis. In this example, if the number of modes of the dual reduced-basis is too small (i.e. smaller than 24), the error estimator does not account for the improvement of u HR when the reduced-basis (ψ k ) N ψ k=1 has more than 3 modes. The convenient choice of N σ ψ is set-up by considering its influence on the smallest value of η for large values of N ψ . Here, N σ ψ = 29 makes this influence negligible for all values of N ψ . When N σ ψ has been fixed, we can estimate c η according to Equation (66). Figure 5 shows the curve related to the HRAE and the error estimator c η η obtained for the sampling point μ = μ 1 , when varying the number of displacement modes for the maximization problem (66).

For the proposed example, inside the plastic zone, the stress field through the thickness of the beam is very sensitive to modifications of the parameters h and C. Here, the validity domain of the hyper-reduction approximations is defined by D = 0.05 2 N t n=1 (u n exp ) 2 . It is shown in Figure 6 step. The speed-up obtained with this hyper-reduced model, including the error estimation, is 6. Without error estimation the speed-up is 8. Hence, the cost of error estimation is quite reasonable and the computational-complexity reduction is preserved. Here, the FE model is very simple. The bigger the FE model, the larger we expect the speed-up to be. The estimation of the validity domain is shown on the right of Figure 6. It is quite restrictive. This means that most of the points in D V are in the reference validity domain D V .

Figure 6 Error prediction on outputs. On the left, in blue, the reference validity domain D V over the parameter space D, on the right the approximate one D V . The white star is located at the sampling point μ 1 . The black dashed line is the limit of the reference validity domain defined by an error on the outputs lower than 5%.

Here, for small values of the parameter h, the prediction of the validity domain is too conservative. Such a situation can appear when the average influence of a parameter on the solution is larger than its influence on the output, because the error estimator is linked to the HRAE in an average sense.

We have tried the proposed numerical approach in case of N ψ = N c ψ (here N c ψ = 6), by using the same RID. The estimated validity domain and the boundary of the reference validity domain for the related HR model are shown in Figure 7. This modification has significant influence on the numerical results. Once again, all the points in D V are in the reference validity domain.

We have performed additional numerical simulations by substituting a Dirichlet boundary condition to the pointwise load. In this case, the time evolution of the load is:

u n c = -U max sin π n N t y, U max = 1.8 mm, N t = 5 (74)
Figure 8 shows the estimation of the validity domain of the hyper-reduced model compared to the boundary of D V , for this case of a non-zero Dirichlet boundary condition. Both reference validity domain and its estimation are different compared to the previous Figure 7 Error prediction on outputs for N ψ = 6. In blue, the approximate validity domain D V . The white star is located at the sampling point μ 1 . The black dashed line is the limit of the reference validity domain defined by an error on the outputs lower than 5%. ones because the reduced bases are different. They have a better accuracy with respect to N ψ = 7 and N σ ψ = 29. We can notice that the higher h and C, the smaller the plastic strains in the beam. Since the hyper-reduced model is accurate for low plastic strains, the validity domain covers the high values of h and C.

Conclusions

We propose an a posteriori estimator of hyper-reduction errors that aims to evaluate if the simulation outputs predicted by hyper-reduced models are convenient approximations of the outputs that the finite element simulation would have predicted. By choosing a tailored norm of the approximation error on displacement, we show how the proposed error estimator is related to the HR approximation error. This error estimator receives the contribution of a coupling term between the error at time step n and the error committed before. We show that this term prevents rigorous error certification. But in practice, numerical examples show a restrictive estimation of the validity domain of the HR model. By restrictive, we mean that this domain is inside the reference validity domain computed by using both the FE solution and the HR solution. The speed-up achieved by the HR model including the estimation of the constitutive relation error is 6 on the proposed numerical elasto-viscoplastic examples. The estimation of the validity domain requires the computation of a constant similar to a Lipschitz constant. This constant is estimated without any additional simulation of FE solution, but the FE solution at the sampling point μ = μ 1 . As shown in is this paper, the numerical implementation of the proposed error estimator is very simple. It does not require the computation of the incremental potential or of its dual. This paper is a first attempt for a posteriori error estimation of hyper-reduction approximations. More numerical experiments are in preparation for the assessment of the proposed numerical approach.

Figure 1

 1 Figure 1 Tailored norm. Schematic view of the physically-based norm of the approximation error, related to the convexity of the incremental potential. The ordinates of points on the dashed line are given by z (ε ) = Nt

Figure 2

 2 Figure 2Meshes. On the top, we show the boundary condition and the mesh of the FE model. On the bottom, we show the RID involving 31 elements of the original mesh.

  71) δ(x = 0) = 0 (72) where E = 80000. MPa, ν = 0.3, m = 8., R o = 0.1 MPa. These material coefficients are reasonably good enough to fit experiments performed at the Centre des Matériaux. The parameters of the mechanical model are the material coefficients h and C related to viscoplasticity and kinematic hardening respectively. The range of variation for μ = (h, C) is [ 15. MPa, 30. MPa] ×[ 3000. MPa, 40000. MPa]. The pointwise load follows a sinusoidal time variation:

  σ ψ = 48 admissible stress modes. We denote by (ψ σ k ) N σ ψ

Figure 3

 3 Figure 3 Empirical modes. Some POD modes: on the left for displacements and on the right for admissible stresses.

4

  . The hyper-reduced model is set-up by choosing HR = D 25. 10 2 . Then, N ψ = 7 and the related RID Z , shown in Figure 2, involves 31 elements over a total of 120 elements. The average number of Newton-Raphson iterations is 4 per time Convergence of the error estimator. The constitutive relation error depends on both displacement approximation and admissible stress approximation. Here, μ = μ 1 .

2 Figure 5

 25 Figure[START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF] Constant fitting. The coefficient c η is fitted upon the simulation outputs obtained for μ = μ 1 such that the error estimator gives an upper bound of the output error, in the case μ = μ 1 .

2 Figure 8

 28 Figure 8Error prediction on outputs in case of non-zero Dirichlet boundary condition. In blue, the approximate validity domain D V . The white star is located at the sampling point μ 1 . The black dashed line is the limit of the reference validity domain defined by an error on the outputs lower than 5%.

Acknowledgements

We would like to thank Pierre Ladevèze, Gérard Coffignal and Yvon Maday for the discussions we had about this work.

Abbreviations

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

The theory represents joint work by DR and LG. Numerical results represents joint work by DR and SJ. All authors read and approved the final manuscript.