Wouter J T Bos

Jean-Pierre Bertoglio

Large-scale bottleneck effect in two-dimensional turbulence

The bottleneck phenomenon in three-dimensional turbulence is generally associated with the dissipation range of the energy spectrum. In the present work, it is shown by using a twopoint closure theory, that in two-dimensional turbulence it is possible to observe a bottleneck at the large scales, due to the effect of friction on the inverse energy cascade. This large-scale bottleneck is directly related to the process of energy condensation, the pile up of energy at wavenumbers corresponding to the domain size. The link between the use of friction and the creation of spacefilling structures is discussed and it is concluded that the careless use of hypo-friction might reduce the inertial range of the energy spectrum.

I. INTRODUCTION

In high Reynolds number three-dimensional turbulence at large wavenumbers k, the compensated energy spectrum k 5/3 E(k) often shows a little positive bump, before it drops down rapidly in the dissipation range. This pre-dissipative bump region was first observed in closure computations [START_REF] André | Influence of helicity on the evolution of isotropic turbulence at high Reynolds number[END_REF], but no physical explanation was given. The bump was also observed in experimental work [START_REF] Mestayer | Local isotropy and anisotropy in a high-Reynolds-number turbulent boundary layer[END_REF][START_REF] Saddoughi | Local isotropy in turbulent boundary layers at high Reynolds number[END_REF] and direct numerical simulations [START_REF] Kerr | Velocity, scalar and transfer spectra in numerical turbulence[END_REF]. A physical explanation of the effect was first given in the work by Herring et al. [START_REF] Herring | A comparative assessment of spectral closures as applied to passive scalar diffusion[END_REF]. They explain that, because the wavenumbers above the Kolmogorov scale are damped by viscosity, some of the nonlocal interactions are damped, which causes the transfer of the energy flux to be less effective, leading to a pile-up of energy. A more detailed investigation of this effect was performed by Falkovich [START_REF] Falkovich | Bottleneck phenomenon in developed turbulence[END_REF], who predicted its wavenumber-dependence, and who named it the bottleneck. In that work it is also inferred that the bottleneck should become more pronounced when hyperviscosity is used. This was illustrated in threedimensional turbulence [START_REF] Borue | Forced three-dimensional homogeneous turbulence with hyperviscosity[END_REF] and two-dimensional turbulence with a direct energy cascade [START_REF] Biskamp | Nonlocal bottleneck effect in two-dimensional turbulence[END_REF] such as magnetohydrodynamic turbulence.

In the present work, we show that in two-dimensional turbulence it is possible to observe a bottleneck at the large scales, due to the effect of friction on the inverse energy cascade. The possibility of large-scale bottlenecks was advanced by Lohse and Mueller-Groeling [START_REF] Lohse | Bottleneck effects in turbulence: scaling phenomena in r versus p space[END_REF]. This range is in three-dimensional turbulence generally affected by a forcing which supplies energy to the flow, complicating the observation of a possible large-scale bottleneck.

An alternative interpretation of the bottleneck in three-dimensional turbulence was recently proposed by Frisch and coworkers [START_REF] Frisch | Hyperviscosity, Galerkin truncation and bottlenecks in turbulence[END_REF]. They relate the bottleneck to statistical equilibrium, and view it as an incomplete thermalization. Indeed, it was stated that if the Laplacian ∇ 2 in the viscous term of the Navier-Stokes equations is replaced by a hyperviscous term ∇ 2α , with the dissipativity α tending to infinity, a Galerkin truncated system is obtained. In Galerkin-truncated 3D inviscid turbulence, energy piles up at the highest wavenumbers, leading to wavemodes in equipartition at the tail of the spectrum. The simultaneous observation of a k -5/3 Kolmogorov inertial range energy spectrum and this equipartition range, proportional to k 2 , was recently observed in direct numerical simulations [START_REF] Cichowlas | Effective dissipation and turbulence in spectrally truncated euler flows[END_REF] (DNS) and closure calculations [START_REF] Bos | Dynamics of spectrally truncated inviscid turbulence[END_REF]. In between these two ranges a pseudodissipation range was observed, caused by nonlocal interactions between the thermal bath and the end of the inertial range.

Since we are interested in the large scales in the present work, we discuss the influence of infrared wavenumber truncation, due to the finite domain-size. For threedimensional turbulence, this truncation influences decaying turbulence and increases the energy decay rate, if the size of the integral lengthscale becomes comparable to the domain-size [START_REF] Touil | The decay of turbulence in a bounded domain[END_REF]. Steady turbulence, forced at a wavenumber k i > k 0 , in which k 0 is the smallest wavenumber of the domain, called the fundamental, will establish an infrared region proportional to k 2 between k 0 and k i , independent of the exact choice of k 0 [START_REF] Lesieur | Turbulence in fluids[END_REF].

In two space-dimensions, the effect of infrared truncation is very different. The absence of vortex stretching induces a double cascade with energy cascading to the large scales and enstrophy going to the small scales. In this case, the large scales are not directly influenced by forcing. The energy spectrum is not necessarily decreasing for k tending to zero, if no energy-sink is present at low k. In an infinite domain, the energy will continue to cascade to smaller and smaller wavenumbers. In the case of a finite domain-size, the energy, cascading backwards, piles up at the fundamental. This condensation was predicted by Kraichnan [START_REF] Kraichnan | Inertial ranges in two-dimensional turbulence[END_REF] and observed numerically [START_REF] Hossain | Long-time states of inverse cascades in the presence of a maximum length scale[END_REF]. The condensation process was studied in detail by Smith and Yakhot [START_REF] Smith | Bose condensation and smallscale structure generation in a random force driven 2d turbulence[END_REF][START_REF] Smith | Finite-size effects in forced two-dimensional turbulence[END_REF], who showed that the statistics of the velocity field remain close to Gaussian until the condensate starts to form. At later times the statistics become non-Gaussian. We note that the condensate state is not the only possible state of finite-size two-dimensional turbulence in the absence of an infrared energy sink. Indeed Tran and Bowman [START_REF] Tran | On the dual cascade in two-dimensional turbulence[END_REF][START_REF] Tran | Robustness of the inverse energy cascade in two-dimensional turbulence[END_REF] showed that if the in-verse energy cascade is weak, i.e., if only a small part of the injected energy cascades to small wavenumbers, the rest being dissipated around the forcing scale, no condensation takes place, but a gradual steepening of the spectrum occurs, when the cascade reaches the smallest wavenumbers, correponding to the box-size. This steepening results in a k -3 range at late times. This scenario was also observed by Chertkov et al. [START_REF] Chertkov | Dynamics of energy condensation in twodimensional turbulence[END_REF]. In the present work we consider the case in which the energy cascade is strong and in which condensation takes place at the fundamental, when no energy sink is present.

II. THE LINK BETWEEN HYPOFRICTION AND ENERGY CONDENSATION

To avoid condensation in two-dimensional turbulence, large-scale friction is often used. Friction acts as an energy sink. The influence of friction, which in principle acts throughout the spectrum, can be concentrated at the smallest wavenumbers by using hypofriction, analogous to hyperviscosity at the highest wavenumbers. The Navier-Stokes equations in this case become:

∂ ∂t + (u • ∇)u = -∇p + ν∇ 2 u + f -µk 2α f (-∇ 2 ) -α u; ∇ • u = 0, (1)
in which u is the velocity vector, p the pressure, ν the viscosity, f an isotropic forcing term localized in spectral space around wavenumber k i , µ the friction coefficient, k f the friction wavenumber, and we shall call α the frictionality, in analogy with the dissipativity in hyperviscous flows [START_REF] Frisch | Hyperviscosity, Galerkin truncation and bottlenecks in turbulence[END_REF].

In the presence of hypofriction, an energy sink is present at low k (e.g. [START_REF] Borue | Inverse energy cascade in stationary twodimensional homogeneous turbulence[END_REF][START_REF] Huang | Anisotropic spectra in two-dimensional turbulence on the surface of a rotating sphere[END_REF]). In this case a steep fall-off of the energy spectrum at low k can be created for large α. If the solution of the hypofrictional case converges to a Galerkin truncated case if α → ∞, in analogy with the work of Frisch et al. [START_REF] Frisch | Hyperviscosity, Galerkin truncation and bottlenecks in turbulence[END_REF], a pile-up of energy can be expected around k f . In the case of linear friction (α = 0), no noticable bottleneck is observed (e.g. [START_REF] Paret | Intermittency in the twodimensional inverse cascade of energy: Experimental observations[END_REF]). The spectral slope d ln(E(k))/d ln(k), is a continuously decreasing function of k for k 0 < k ≪ k i with k 0 the fundamental and k i the wavenumber where energy is introduced by the forcing.

In the present study, using a closure model, we will systematically vary the order of the frictionality from α = 0 to α = 128. We will investigate the appearance of bottlenecks, its influence on the inertial range scaling and the convergence of this bottleneck to an infrared wavenumber truncation.

Some evidence for large-scale bottlenecks can be found in literature [START_REF] Borue | Inverse energy cascade in stationary twodimensional homogeneous turbulence[END_REF][START_REF] Huang | Anisotropic spectra in two-dimensional turbulence on the surface of a rotating sphere[END_REF]. In these studies hypofriction was used with values 8 and 5 respectively. The compensated energy spectra do show a pile up of energy at the small wavenumbers. A link with the creation of coherent structures is made but no systematic study of the influence of α on the large-scale bottleneck was reported in these studies.

Our interpretation, in the light of the present study, is that the bottlenecks in these works are in fact incomplete condensates, in analogy with the interpretation of Frisch et al. [START_REF] Frisch | Hyperviscosity, Galerkin truncation and bottlenecks in turbulence[END_REF], that the bottleneck in 3D turbulence is an incomplete thermalization. The coherent structures that are observed are closely related to the structures, or condensates, that would be observed in a truncated solution.

Even though numerical resources increase rapidly and two-dimensional turbulence is computationally accessible at high Reynolds numbers, a systematic DNS study of the influence of α would be very expensive. Furthermore the small wavenumber range of an energy spectrum converges more slowly to a statistically steady state than the small scales. It was shown that for this kind of studies, twopoint closure is a convenient tool [START_REF] Frisch | Hyperviscosity, Galerkin truncation and bottlenecks in turbulence[END_REF][START_REF] Bos | Dynamics of spectrally truncated inviscid turbulence[END_REF]. For the case of two-dimensional turbulence, it was pointed out [START_REF] Herring | Comparison of direct numerical simulation of two-dimensional turbulence with two-point closure: the effects of intermittency[END_REF] that in the absence of coherent vortices, closure approaches compare well to direct numerical simulation. The inverse cascade is known to be free from large coherent structures and its statistics are close to Gaussian [START_REF] Boffetta | Inverse energy cascade in two-dimensional turbulence: deviations from Gaussian behavior[END_REF]. In direct numerical simulations, if friction is absent, the pile up of energy at the fundamental will eventually give rise to the creation of two counter-rotating vortices of the size of the domain.

In the long-time limit, the statistics become highly non-Gaussian [START_REF] Smith | Bose condensation and smallscale structure generation in a random force driven 2d turbulence[END_REF] so that the use of the model, which are based on an expansion about Gaussianity, will probably not be justified [START_REF] Kraichnan | Is there a statistical mechanics of turbulence?[END_REF][START_REF] Chen | Non-gaussian statistics in isotropic turbulence[END_REF][START_REF] Bertoglio | EDQNM closure: A homogeneous simulation to support it. a quasihomogeneous simulation to disprove it[END_REF]. The time before this happens is supposed to be reasonably well predicted by closure. It is this time-interval, in which pile-up of kinetic energy at small wavenumbers takes place, and which precedes the formation of coherent structures, on which we concentrate in the present work.

III. OBSERVATION OF LARGE-SCALE BOTTLENECKS

We use the Eddy-Damped Quasi Normal Markovian (EDQNM) closure. It was used to study two-dimensional turbulence by Leith [START_REF] Leith | Atmospheric predictability and twodimensional turbulence[END_REF] and Pouquet et al. [START_REF] Pouquet | Evolution of high Reynolds number two-dimensional turbulence[END_REF]. In this closure, the evolution equation for the energy spectrum is solved

∂ t + νk 2 + µ(k f /k) 2α E(k, t) = T N L (k, t) + F (k). (2)
The nonlinear transfer T N L is

T N L (k, t) = 4 π ∆ Θ kpq xy -z + 2z 3 √ 1 -x 2 k 2 pE(p, t)E(q, t) -kp 2 E(q, t)E(k, t) dpdq pq . ( 3
)
The forcing is a time-independent energy input of rate c,

localized at k = k i , F (k) = c δ(k -k i ).
In equation ( 3), ∆ is a band in p, q-space so that the three wave-vectors k, p, q form a triangle. x, y, z are the cosines of the angles opposite to k, p, q in this triangle. The characteristic time Θ kpq is defined as:

Θ kpq = 1 -exp(-(η k + η p + η q ) × t) η k + η p + η q , (4)
in which η is the eddy damping. In the present study we use the following expression for the damping coefficient

η k = λ k 0 s 2 E(s, t)dS + νk 2 + µ(k f /k) 2α . (5)
In expression [START_REF] Herring | A comparative assessment of spectral closures as applied to passive scalar diffusion[END_REF], three contributions to the eddy damping are present. The first two are the standard contributions introduced in the EDQNM closure. For λ we use the value 0.4 [START_REF] Lesieur | Diffusion of a passive scalar in two-dimensional turbulence[END_REF]. The third contribution is an adhoc modification of the damping to take into account the presence of friction at large scale. A rigorous derivation of the modification induced by the friction as done for the rotating case [START_REF] Bellet | Wave turbulence in rapidly rotating flows[END_REF] or using an improved eddy-damping [START_REF] Bos | A single-time two-point closure based on fluid particle displacements[END_REF] is not attempted in the present work. The numerical set-up is similar to the one in [START_REF] Bos | Dynamics of spectrally truncated inviscid turbulence[END_REF]. The wavenumbers are spaced as k p = k 0 r p-1 with k 0 the fundamental and r > 1 a constant which determines the number of wavenumbers per decade. The spatial resolution for the computations reported here is 40 wavenumbers per decade which corresponds to a resolution of F ≈ 12 wavenumbers per octave. Simulations are reported for α = 0, 1, 2, 4, 8, 32, 128 as well as a friction-free simulation truncated at k 0 = k f . For all runs but this last one, a steady state was obtained, for which the energy spectrum is shown in Figure 1. For the friction-free (truncated) case, a steady state will not be obtained until the energy-pile up at the fundamental is balanced by the viscous dissipation at this wavenumber, which is very small in the present case as the Reynolds number is large. The fundamental would then accumulate an amount of energy which is enormous compared to the cases with small frictionality. For a better comparison between the different cases, in the following a spectrum is shown for which the energy pile-up at the fundamental is comparable to the pile-up for α = 128. In Figure 2 we observe that for α > 0 a bottleneck is present around k f ≈ 3. The size of this bottleneck increases with α. For values of 1 ≤ α < 8, the bottleneck has the form of a bump. For α ≥ 8, the bottleneck changes its shape towards a peak and for the highest values of α used in the present work, the bottleneck collapses with the energy spectrum of the frictionless truncated case. The k -5/3 inertial range is influenced by this bottleneck, not only on the bottleneck itself, where the spectrum is steeper than k -5/3 , but also in the region with wavenumbers larger than the bottleneck, and extending for about one decade, the spectrum deviates and is slightly shallower than k -5/3 . The friction-free truncated case, and the high-α friction cases, show much similarity (at small k) with the inviscid truncated 3D case (at high k) studied in [START_REF] Bos | Dynamics of spectrally truncated inviscid turbulence[END_REF]. A depletion of the compensated energy spectrum is present in the vicinity of the peak (for k > k f in the 2D case, and for k < k th in the 3D case, where k th is the wavenumber marking the beginning of the equipartition range). The effect can, in both cases, be explained by nonlocal interactions, as was done in the three-dimensional case [START_REF] Bos | Dynamics of spectrally truncated inviscid turbulence[END_REF].

It is interesting to investigate whether this depleted region in the compensated energy spectrum corresponds to the emergence of a new scaling regime. Dimensional analysis based on the assumption of a constant flux of energy and a sweeping timescale related to the advection of the small eddies by the velocity associated with the kinetic energy contained in the condensate, one can deduce a scaling proportional to k -3/2 . It is observed in the inset of Figure 2 that a small wavenumber range shows a spectral exponent close to this value (1.53). This range is however too short to make conclusive statements.

It was shown by Smith and Yakhot [START_REF] Smith | Bose condensation and smallscale structure generation in a random force driven 2d turbulence[END_REF], that when the condensate becomes large enough, the statistics of the velocity field become strongly non-Gaussian. One could then ask, if the present closure approach, based on an assumption of maximum randomness [START_REF] Kraichnan | The structure of isotropic turbulence at very high Reynolds numbers[END_REF], is still valid. To estimate the extend to which the closure gives results for the energy spectrum close to results of DNS, we compare our results with DNS results from literature [START_REF] Hossain | Long-time states of inverse cascades in the presence of a maximum length scale[END_REF]. In that study the coexistence of a large peak at the fundamental and a k -5/3 range for larger wave-numbers is observed. The value of the energy accumulated at the fundamental is one order of magnitude above its inertial range value. In Figure 2 the peak value is also one order of magnitude above its inertial range value. So at least qualitatively, the coexistence of a large peak and the k -5/3 inertial range is observed in DNS as well as in closure. However this validation remains qualitative and the present results could only be fully supported by a detailed and quantitative comparison with DNS. Such a comparison is beyond the scope of the present paper.

When the bottleneck becomes large enough, a small kink is observed in the spectrum around k = 50. This comes from the discretization. Nonlocal interactions are not resolved if the ratio of the smallest wavenumber to the middle wavenumber of an interacting triad is smaller than [START_REF] Pouquet | Evolution of high Reynolds number two-dimensional turbulence[END_REF] a = 2 1/F -1. For elongated triads, the middle wavenumber and the longest wavenumber are approximately of the same size. The largest wavenumber k directly interacting with k f can then be shown to be k ≈ 50 for F = 12. This clearly illustrates the nonlocal origin of the pseudo-friction range, created by the bottleneck. In principle, the effect of the truncation can be removed by refining the discretization. To allow nonlocal interactions over the full 4 decades of the simulation, the number of necessary gridpoints can be estimated to be roughly 700 per octave which would make the simulations prohibitively expensive. Another solution would be to explicitly compute the first-order contributions of the nonlocal interactions [START_REF] Pouquet | Evolution of high Reynolds number two-dimensional turbulence[END_REF] by expanding the nonlinear transfer with respect to the parameter a, which is small for very nonlocal interactions. This procedure has the disadvantage not to conserve energy or to mix-up the different orders of the expansion. The improvement expected to be minimal with respect to the presented results, these solutions will not be attempted here.

IV. CONCLUSION

The present results show that the bottleneck effect is not only observable at the small scales, but it can also be created, be it artificially, at the large scales. It shows that by increasing α, the energy spectrum tends to the condensed solution predicted by Kraichnan [15]. This constitutes a logical extension of the ideas of Frisch et al. [START_REF] Frisch | Hyperviscosity, Galerkin truncation and bottlenecks in turbulence[END_REF] to the inverse cascade in two-dimensional turbulence. It also shows that pile-up of energy at small wavenumbers influences the inertial range scaling by nonlocal interactions. The present work thereby confirms the ideas by Sukoriansky et al. [START_REF] Sukoriansky | Large scale drag representation in simulations of twodimensional turbulence[END_REF][START_REF] Sukoriansky | Subgrid-and supergrid-scale parameterization of turbulence in quasitwo-dimensional barotropic flows and the phenomenon of negative viscosity[END_REF], that large-scale friction can influence all scales of two-dimensional turbulence, since the pile-up of energy at the small scales might trigger the creation of coherent structures. These coherent structures can induce a k -3 scaling [START_REF] Chertkov | Dynamics of energy condensation in twodimensional turbulence[END_REF][START_REF] Borue | Inverse energy cascade in stationary twodimensional homogeneous turbulence[END_REF], which camouflates the inertial range. The use of hypofriction, intended to enlarge the inertial range, can thus have the opposite effect.

A detailed DNS study would be necessary to confirm the effects observed in the present study. Even at low or moderate resolution, trends could probably be observed and answer the question whether the predictions of the present closure, based on an expansion about Gaussianity, are valid for the flow here considered in which coherent large scale structures are present. In the authors' opinion, in this type of flow, as is the case for many types of turbulent fields, the simultaneous use of Direct Numerical Simulations at low and moderate Reynolds numbers and of closures at high Reynolds numbers will probably constitute the best strategy to gain understanding of the physics of turbulence dynamics. Such a study is however left for future work. Pauls, Samriddhi Sankar Ray, Jian-Zhou Zhu and two anonymous referees.

FIG. 1 :

 1 FIG.1: Steady-state energy spectra with frictionality ranging from 0 to 128 and a friction free case truncated at the low wavenumbers at k = k f .

FIG. 2 :

 2 FIG.2:The same spectra as in Figure1compensated by k 5/3 and, in the inset, compensated by k 1.53 .

Acknowledgments

We acknowledge interaction and discussion with Claude Cambon, Uriel Frisch, Fabien Godeferd, Walter


R Rubinstein

W J T Bos

On the unsteady behavior of turbulence models

Periodically forced turbulence is used as a test-case to evaluate the predictions of two-equation and multiple-scale turbulence models in unsteady flows. The limitations of the two-equation model are shown to originate in the basic assumption of spectral equilibrium. A multiple-scale model based on a picture of stepwise energy cascade overcomes some of these limitations, but the absence of nonlocal interactions proves to lead to poor predictions of the time variation of the dissipation rate. A new multiple-scale model that includes nonlocal interactions is proposed and shown to reproduce the main features of the frequency response correctly.

A basic premise of one point closures such as the k -ǫ model is the hypothesis of "spectral equilibrium," which justifies two distinct roles of the dissipation rate ǫ: on the one hand, it appears in the energy balance, defined as a correlation of velocity gradients, hence a small-scale quantity; on the other hand, it is used phenomenologically to describe large-scale transport properties. The most basic formulation of the latter is Kolmogorov's hypothesis ǫ ∝ k 3/2 /L where k is the turbulent kinetic energy and where L is a length scale characteristic of the largest scales of motion; equivalent formulations include ν t ∝ k 2 /ǫ, the formula for turbulent viscosity, or τ ∝ k/ǫ, the formula for the turbulent time-scale. The Kolmogorov theory, or more general assumptions of selfsimilarity of all scales of motion, justify all of these proportionalities [START_REF] Oberlack | Similarity in non-rotating and rotating turbulent pipe flows[END_REF], although the constants of proportionality need not coincide in all self-similar flows [START_REF] Bos | Spectral imbalance and the normalized dissipation rate of turbulence[END_REF].

But turbulence models are not needed to describe selfsimilar flows, which merely serve as calibration cases; models are needed to describe departure from selfsimilarity, when spectral equilibrium becomes a strong constraint on turbulence evolution. In a study of a flow in which turbulence evolves from a steady state to a selfsimilar time-dependent state [START_REF] Rubinstein | Timedependent isotropic turbulence[END_REF], the implications of the departure from spectral equilibrium were investigated: this departure was connected to transient failure of the Tennekes-Lumley balance [START_REF] Tennekes | A first course in turbulence[END_REF] and to the consequent relevance of small scale dynamics for the large scales. The breakdown of spectral equilibrium has also been observed in engineering flows including turbulent diffusers and wakes. [START_REF] Touil | A spectral closure applied to anisotropic inhomogeneous turbulence[END_REF] The limitations of the spectral equilibrium hypothesis are very well known in the modeling literature and have led to proposals for multiple-scale models [START_REF] Schiestel | Multiple-time-scale modeling of turbulent flows in one-point closures[END_REF] that more realistically address the complexity of the nonlinear interactions in turbulent flows. This letter reports on some investigations of multiple-scale models applied to an especially simple and attractive test case for transient turbulence: periodically forced turbulence [START_REF] Lohse | Periodically kicked turbulence[END_REF][START_REF] Cadot | Observation of resonances in modulated turbulence[END_REF][START_REF] Kuczaj | Turbulence modification by periodically modulated scale-dependent forcing[END_REF]. This problem arises when isotropic incompressible turbulence, maintained in a steady state by a large-scale isotropic forcing with total amplitude p, is subjected to a small time-dependent periodic perturbation with amplitude p: p(t) = p + p cos(ωt), such that the ratio p/p ≪ 1 and such that the forcing length scale does not depend on time. The phase-averaged kinetic energy k can then be decomposed into a mean k and a periodic part k(ω) cos(ωt + φ k (ω)), with φ k the phaseshift between the forcing and the kinetic energy. Similarly, the viscous dissipation rate can be written as ǫ = ǭ + ǫ(ω) cos(ωt + φ ǫ (ω)); k and ǭ are related to the time-independent forcing length scale L by L ∝ k3/2 /ǭ. The periodic parts of k and ǫ are sinusoidal, like the forcing, because p/p ≪ 1. The functions k(ω), ǫ(ω), φ k (ω), and φ ǫ (ω), which characterize the linearized response of steady-state turbulence to periodic perturbation of the forcing, can be called the linear response functions.

We will use periodically forced turbulence as a test case to evaluate the ability of multiple-scale models to predict the dynamics of time-dependent turbulence. It will first be shown that the unsteady predictions of multiple-scale models are significantly better than the predictions of a two-equation model. However, an elementary multiplescale model based on the heuristic picture of stepwise energy cascade is found to have limitations in predicting the unsteady dissipation rate. A multiple-scale model that includes the possibility of nonlocal interactions is proposed; it is shown that this model can capture some fine features of the unsteady energy dissipation.

The linear response functions were determined in recent work [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] using the EDQNM closure theory, which was shown to compare very well to available low Reynolds number experimental, and DNS data. Comparison with high Reynolds number data for linear response functions would be desirable, but such data is not yet available. Briefly summarizing the major conclusions, the two-equation model is satisfactory both in the static limit ω → 0, in which the phase shifts φ k , φ ǫ vanish, and in the frozen limit ω → ∞, in which k ∼ ω -1 and φ k ∼ π/2. However, this agreement is trivial; only the results at intermediate frequencies provide a real test of the model. The two-equation model reproduces the function k(ω) reasonably well, but the transition from the static limit φ k ≈ 0 to the frozen limit φ k ≈ π/2 occurs over a frequency range that is much too wide. The amplitude ǫ(ω) is not satisfactory: a range in which ǫ ∼ ω -3 at high Reynolds numbers is absent. We note that since observations of the modulated dissipation rate are very difficult, and relevant high Reynolds number data is not yet available, the EDQNM results for this quantity remain theoretical predictions; they are nevertheless supported by arguments [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] based on the well-established role of distant interactions in turbulence. Finally, the two-equation model also makes the incorrect prediction that φ k = φ ǫ regardless of ω.

We investigate whether these predictions can be improved using multiple-scale modeling following the ideas of Schiestel [START_REF] Schiestel | Multiple-time-scale modeling of turbulent flows in one-point closures[END_REF]. Multiple-scale models can be considered numerical methods for spectral closures, with significant modifications designed to permit reasonable accuracy at a very low order of discretization. Thus, whereas a numerical implementation of a spectral closure would solve for the energy spectrum at perhaps hundreds of discrete wavenumbers κ i , in Schiestel's formulation, the energy spectrum is divided into a relatively small number of wavenumber shells κ i-1 ≤ κ ≤ κ i ; for each shell, equations are written for two scalar descriptors: the fluctuation energy contained in the shell and the net energy flux into it. To enhance the accuracy possible with a relatively small number of shells, Schiestel allowed the partition wavenumbers κ i to be functions of time. A schematic picture of the resulting discretized energy cascade is given in Figure 1, following Lumley [START_REF] Lumley | Some comments on turbulence[END_REF].

The starting point for the analytical formulation is the Lin equation governing the energy spectrum E(κ, t),

∂ ∂t + νκ 2 E(κ, t) = P (κ, t) - ∂F (κ, t) ∂κ . (1)
In this equation ν is the viscosity, P (κ, t) the forcing term, and F (κ, t) is the energy flux across wavenumber κ. Schiestel applied the Kovaznay model [START_REF] Monin | Statistical Fluid Mechanics[END_REF]

F (κ, t) = Cκ 5/2 E(κ, t) 3/2 , (2)
with C a parameter which determines the Kolmogorov constant. This model represents a stepwise cascade of energy in spectral space from small to large κ and reproduces a κ -5/3 inertial range. Figure 1 would correspond to a stepwise cascade if the dashed lines were absent.

To obtain a multiple-scale model, write the equation for the time derivative of the spectral energy flux F (κ, t),

Ḟ (κ, t) = 3 2 F (κ, t) E(κ, t) Ė(κ, t). (3)
From the viewpoint of Schiestel's analysis, we have assumed that the partition wavenumbers are constant in time: this assumption seems appropriate for this problem, in which the forcing wavenumber is fixed. The production scales are assumed to be confined to low κ, and the dissipation to high κ (a more general description including finite-Reynolds number effects or broad-band forcing will not be attempted here). For κ in the inertial range, (1) becomes

Ė(κ, t) = - ∂F (κ, t) ∂κ . (4)
We combine this with ( 2) and ( 3) to obtain

Ḟ (κ, t) = - 3 2
F (κ, t) E(κ, t) ∂F (κ, t) ∂κ . (5)
A discrete model is obtained by splitting the spectral domain into n shells as illustrated in Figure 1. The energy in shell i is e i ≈ E(κ i )∆κ i . The spectral flux f i ≈ F (κ i ) and the time derivative is ∂F (κ i , t)/∂k ≈ (f i -f i-1 )/∆κ i , so that we obtain ḟi = -3 2

f i e i (f i -f i-1 ), 1 ≤ i ≤ n. (6)
Integrating ( 4) over each shell gives the partial energy balance equations

ėi = -(f i -f i-1 ), 1 ≤ i ≤ n. (7)
in which all f i and e i are functions only of time. The oscillating production term p, assumed to act at the small wavenumbers, is identified with f 0 , the flux entering shell 1. Furthermore, the high Reynolds numbers case is considered in which we assume that the viscous dissipation takes place at the last wavenumber shell: ǫ = f n . This assumption makes it unnecessary to introduce partial dissipation rates for each shell and corresponding equations of motion. A special feature of the Kovaznay model is that the partition wavenumbers do not appear in the model. By choosing n = 1, one obtains a two-equation model; the equation for f 1 becomes the dissipation rate equation ǫ = (3/2)(ǫ/k)(P -ǫ). Note that the two model constants, generally called C ǫ1 and C ǫ2 in the literature, are equal, which allows the study of statistically stationary isotropic turbulence. Consistency with homogeneous shear flow, or with any problem in which the forcing length scale increases as a power law or exponential in time [START_REF] Rubinstein | Self-similar turbulence evolution and the dissipation rate transport equation[END_REF], requires C ǫ1 < C ǫ2 .

Choosing n ≥ 2 should improve the predictions by introducing the possibility of spectral imbalance, a necessary requirement if the same model is to be applied to both forced and decaying turbulence [START_REF] Bos | Spectral imbalance and the normalized dissipation rate of turbulence[END_REF]; imbalance is possible because the partial fluxes f i with 1 ≤ i ≤ n -1 and the dissipation ǫ = f n are independent.

The unsteady behavior of this model is now assessed as follows. Beginning with a steady state with shell energies ēi such that p = fi = ǭ, the periodic perturbation p cos(ωt) is added to the forcing. Describing the periodic response in terms of complex amplitudes k(ω) and ǫ(ω), so that k(ω) = | k(ω)| and tan φ k (ω) = ℑ k(ω)/ℜ k(ω) with the obvious analogs for ǫ, the equations for the periodic part of the partial kinetic energies and local fluxes become

iωê i = -( fi -fi-1 ), iω fi = - 3 2 fi ēi ( fi -fi-1 ). ( 8
)
The resulting linear system is easily solved analytically for êi (ω) and fi (ω) in terms of the periodic forcing perturbation p and the parameters ēi and fi . Figure 2 compares the results for the amplitude k(ω) and phase shift φ k (ω) for models with n = 1 to n = 7 wavenumber shells. Also shown are results obtained in EDQNM computations [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] at Reynolds number R λ = 1000, with R λ ≈ 15R L 1/2 and R L = (2 k/3) 1/2 L/ν. It has been shown [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] that at low ω, k(ω), should tend to a plateau and that at large ω, k(ω) follows a powerlaw proportional to ω -1 . Confirming the conclusion of Bos et al. [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF], k(ω) is in reasonable agreement with EDQNM even for the two-equation model n = 1, although this agreement improves significantly as the number of wavenumber partitions n increases.

The error in the two-equation model prediction of k(ω) is a too gradual transition from the plateau to the ω -1 region. This defect appears more prominently in the phase φ k (ω): although all models give the correct static and frozen limits, the two-equation model transitions much too gradually, and only the models with n ≥ 2 are in close agreement with EDQNM. The relatively rapid transition in the energy phase shift therefore appears as a typical multiple-scale effect: evidently, in this problem, the small scales are not simply 'slaved' to the large scales through a constant dissipation rate as is assumed in a two-equation model; instead, they are dynamically independent and have a strong effect on what is apparently a purely largescale property.

It has also been demonstrated [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] that the response function ǫ(ω) follows at high ω a power-law proportional to ω -3 up to the Kolmogorov frequency ω η ∼ ǭ/ν. For ω > ω η , ǫ(ω) becomes proportional to ω -1 . These results are shown for comparison in the graph on the left side of Figure 3 (the phase shift φ ǫ proves difficult to compute with any confidence because of the extremely small amplitudes involved, therefore comparisons are omitted).

The agreement of ǫ(ω) given by the multiple-scale model Eq. ( 8) with the EDQNM results is very good down to values ǫ/p = 10 -3 for n > 3 (note that k(ω) and ǫ(ω) are proportional to p, and p is chosen unity without loss of generality). However for smaller values of ǫ/p, the discrete model starts to diverge from the EDQNM results, especially for large n. It is very easily shown that for this model, the leading order contributions at high ω are proportional to ω -n : indeed, recursive solution of the equations for fi in [START_REF] Cadot | Observation of resonances in modulated turbulence[END_REF] show that fi ∼ ω -i p. Thus, only if n = 3 can we obtain ǫ(ω) ∼ ω -3 , but this is the result of a coincidence, which disappears if the number of partition wavenumbers is increased. This difficulty reflects a limitation of the multiple-scale model: Eq. ( 5) implies a linear first-order partial differential equation for F in which disturbances in F propagate along characteristics; this property is probably significantly compromised by a finite dimensional approximate model. It has been shown [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] that nonlocal interactions are responsible for the ω -3 range. Nonlocal interactions are represented in Figure 1 by dashed lines. Such interactions do not occur in the model Eqs. ( 6) and ( 7), because quantities in shell i depend only on its nearest neighbor, shell i -1. The absence of nonlocal interactions will be even more significant in problems in which the role of nonlocality is greater, as in the Batchelor regime of the passive scalar [START_REF] Batchelor | Small-scale variation of convected quantities in turbulent fluid[END_REF] and in MHD [START_REF] Pouquet | Strong MHD helical turbulence and nonlinear dynamo effect[END_REF].

To address this problem, a new multiple-scale model will now be derived including the effect of nonlocal interactions. We start from a simple spectral model containing nonlocal interactions, due to Ellison [START_REF] Monin | Statistical Fluid Mechanics[END_REF],

F (κ, t) = C κ 0 κ 2 E(κ, t)dκ 1/2 κE(κ, t). (9)
Applying the same procedure as to the Kovaznay model, the partial energy balance Eq. ( 7) is unchanged, but Eq. ( 6) is replaced by ḟi = -

f i e i (f i -f i-1 ) - f i 2 i p=1 κ 2 p (f p -f p-1 ) i p=1 κ 2 p e p . (10)
This expression contains the wavenumbers κ i because nonlocal interactions, which depend on the spacing between the wave-number partitions, have been retained. The specification of the κ i becomes part of the model: we will use a logarithmic discretization with κ i = κ 1 r i-1 in which r is a model parameter which determines the logarithmic gridsize. Using this discretization, the ratio κ n /κ 1 = r n-1 , so that a large range of scales can be considered by increasing r. Note that if a linear discretization is used, κ n /κ 1 = (n∆κ/1∆κ) = n so that the number of partitions for high Reynolds numbers becomes prohibitively large. Using the logarithmic discretization, the model for periodic forcing becomes

ḟi = - fi ēi ( fi -fi-1 ) - fi 2 i p=1 r 2(p-1) ( fp -fp-1 ) i p=1 r 2(p-1) ēp , (11)
with the same partial energy equations as in Eq. ( 8). Again, when n = 1 the model reduces to a two-equation model. For n > 1 the model differs from the previous model through the interaction term which couples wavenumber shell i with all wavenumber shells p = 1, • • • , i.

The response functions depend on both n and r, which will be chosen through a compromise between computational cost and precision.

The results for ǫ(ω) obtained from Eq. ( 11) with r = 3 and 1 ≤ n ≤ 7 wavenumber partitions are shown in the graph on the right side of Figure 3. In the same figure, EDQNM results [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF] at R λ = 1000 are shown. The agreement with EDQNM is very good down to values ǫ/p = 10 -3 for 3 < n < 7, and for n = 7, agreement is good down to ǫ/p = 10 -6 . We conclude that the model including nonlocal interactions Eq. ( 11) makes better predictions of ǫ than the model Eq. ( 8), in which nonlocal interactions are absent. The predictions of Eq. ( 11) (not shown) for k(ω) and φ k (ω) very nearly coincide with the results obtained using Eq. (8).

To conclude, we have found that in the problem of periodically forced turbulence, a two-equation model only gives satisfactory predictions for k(ω) and φ k (ω) at asymptotically high and low frequencies. The predictions at intermediate frequencies are improved by using a multiple-scale model based on the heuristic picture of stepwise energy cascade, the Kovaznay model. In particular, this model correctly predicts the rapid jump of the phase shift φ k (ω) between the static and frozen limits. The multiple-scale model based on the Ellison closure includes nonlocal effects, and leads to better agreement with EDQNM, including the high Reynolds number ω -3 scaling range [START_REF] Bos | Small scale response and modeling of periodically forced turbulence[END_REF]. Both models give practically indistinguishable results for the amplitude and phase of the oscillating kinetic energy, which is not strongly influenced by nonlocal interactions. Our approach suggests how one might construct reduced order models for phenomena dominated by significant nonlocal interactions, like the Batchelor range of the passive scalar and some cases of MHD.

We would like to acknowledge the interesting and thoughtful comments of the referees, which led to significant modifications of the paper.

FIG. 1 :

 1 FIG. 1: Schematic view of a discretized energy cascade. The solid arrows indicate the energy fluxes between neighboring wavenumber shells. The dashed lines indicate the nonlocal fluxes between spectrally remote shells. Note that a similar picture can be found in a paper by Lumley [11].

FIG. 2 :

 2 FIG.2: Linear response functions k(ω) (left) and φ k (ω) (right) for the multiple-scale model based on the Kovaznay spectral closure[START_REF] Schiestel | Multiple-time-scale modeling of turbulent flows in one-point closures[END_REF]. The 7 curves in each graph correspond to models with n = 1, 2, .., 7 spectral shells and n increases in the direction of the arrow. The theoretical amplitude prediction ω -1 is indicated by a straight line. Also shown are the results of EDQNM simulations (symbols). The frequency ω is normalized by the large scale frequency ǭ/ k.

FIG. 3 :

 3 FIG. 3: Frequency response function ǫ(ω) for the model based on the Kovaznay closure (8) (left) and for the model based on the Ellison closure (11) (right). The 7 curves in each graph correspond to n = 1, 2, .., 7 and n increases in the direction of the arrow. The theoretical prediction ω -3 is indicated by the straight lines. Also shown are the results of EDQNM simulations (symbols).


Diane Cinca

Study of the top quark pairs production collected by the ATLAS detector at the LHC collider

Interest of Top quark physics

Because of its high mass the top quark is an important milestone to understand Particle Physics. The measurement of its parameters will allow to test the precision of the standard model. The precise measurement if its mass coupled to the one of the W boson will constrain the Higgs boson mass.

Being the heaviest particle discovered it is also a privileged partner for New Physics particles and for the same reason the main background for their search. The precise reconstruction of the Top quark is then a key point for these measurements. In our study we focus on the jet selection for the hadronic top quark.

2 Jet selection in t t → eνb bq q and t t → µνb bq q

At the LHC top quarks will be predominantly produced in pairs due to quarks or gluons fusion. The most promising desintegration channel for our study is the semileptonic one which has a branching ratio about 30%. With only one neutrino, one lepton and four jets in the final state, it is a good compromise between signal efficiency and background level. We focus on the hadronic top quark reconstruction. Giving three jets in the final state, the issue is to select the jets trully coming from the hadronic top with the best efficiency.

The procedure 1 is to select the jets by topological proximity a . The two jets emitted by the W boson are chosen as the two closest in ∆R. When the W boson is reconstructed, the closest jet in ∆R is selected as being the one coming from the b quark. We can then reconstruct the top quark with a good efficiency. Our goal is to improve this method by using multivariate analysis.

Multivariate analysis

On the contrary of the geometrical method which uses only one variable, we want to use many variables combined in multivariate analysis in order to enhance the purity, rate of well reconstructed events against bad, and efficiency, rate of well reconstructed events against events which can be reconstructed, of our sample. Variables implemented are angular parameters, reconstructed W boson mass and jets energies.

Conclusion and perspectives

On figure 1 we use multivariate analysis for the jet selection and see that we improve the efficiency greatly. This method allows us to gain 35% of well reconstructed top quarks compared to geometrical analysis giving us a powerful tool for top quark mass reconstruction. We expect that this method will lead to significant improvements on the top mass resolution.

  Study of the top quark pairs production collected by the ATLAS detector at the LHC collider Diane Cinca Laboratoire de Physique Corpusculaire -Univ. B. Pascal -CNRS/IN2P3 24 Av. des Landais, F-63177 Aubière Cedex -France In nominal conditions the LHC will collide proton beams of 7 TeV producing one top quark each second. I introduce the interest of the top quark physics and an improved procedure for the mass reconstruction in the semileptonic channel developed on Monte Carlo simulations. I show the interest of multivariate analysis for the jet selection.

Figure 1 :

 1 Figure 1: Purity of the sample versus the efficiency of the signal. Black squares show the results for geometrical jet selection. Lines represent the results for purification of the angular jet selected sample. Stars represent the results for MVA jet selection.


D Ryckelynck

D Missoum Benziane

Multi-level a priori hyper reduction of mechanical models involving internal variables

Keywords: APHR method, POD, parallel computing, FETI, domain decomposition, viscoplasticity

This paper concerns the adaptation of reduced-order models during simulations of series of elastoviscoplastic problems. In continuation with previous works, this paper aimed at extending the A Priori Hyper Reduction method (APHR method) for nonlinear thermal problems to nonlinear mechanical problems involving internal variables. This method is an a priori approach because full incremental responses of detailed models are not forecasted in order to build reduced-order models. The recent extension of the Hyper Reduction method to reduction of mechanical models involving internal variables makes possible the reduction of degrees of freedom and the reduction of integration points. A multi-level formulation is introduced to focus on the capability of the method to perform efficient parallel computations to adapt reduced-order models.

Introduction

Elastoviscoplastic models are widely used for fatigue life prediction or crack growth of metallic components [START_REF] Tanguy | Ductile to brittle transition of an A508 steel characterized by Charpy impact test Part II: modeling of the Charpy transition curve[END_REF][START_REF] Nicouleau | Structural calculation and lifetime-prediction in thermomechanical fatigue of engine components[END_REF][START_REF] Mazière | Overspeed burst of elastoviscoplastic rotatig disks -Part I: Analytical and nmerical stability analyses[END_REF]. Usually, the constitutive laws are described in the framework of the irreversible thermodynamic processes. The strain history is taken into account using internal variables. These variables are the lump sum of the history of material changes. This approach has its roots in the works by Biot [4] ; Ziegler [START_REF] Ziegler | Some extremum principles in irreversible thermodynamics with applications to continuum mechanics[END_REF]; Germain [START_REF] Germain | Cours de mecanique des milieux continus[END_REF] or Halpen and Nguyen [START_REF] Halpen | Sur les materiaux standards generalises[END_REF] and has proven its ability to cover a broad spectrum of models in viscoelasticity, viscoplasticity, plasticity and also continuum damage mechanics. Examples of such constitutive laws can be found in [START_REF] Lemaitre | Mecanique des materiaux solides[END_REF]. To simplify the presentation of the multi-level APHR method, we consider small-displacement small-strain mechanical problems having unique solutions. According to the framework of the irreversible thermodynamic processes, a constitutive law can be defined by a choice of: internal variables z, a free energy w(ε, z) and a pseudo potential of dissipation ϕ * [START_REF] Germain | Continuum Thermodynamics[END_REF] in case of standard formulations. Some conjugated variables Z are associated to the internal variables z using a definition of the dissipation. It turns out that the problems related to such mechanical models are nonlinear and time dependent. The Finite Element (FE) method [START_REF] Zienkiewicz | Finite Element Method[END_REF], implicit time integration scheme and the Newton Raphson algorithm are usually used for approximately solving this kind of problem.

The development of large FE models increases the need of low order models created by model reduction methods. The availability of Reduced Order Models (ROMs) can greatly facilitate the solution of series of mechanical problems appearing in optimization problems for instance. The formulation of the reduced-equations differs from the formulation of the detailed equations in the choice of the functional space related to the primal state variables. Most of the time, the state variables are considered as a linear combination of known fields. In case of an a posteriori approach, these fields are obtained by solving detailed preliminary problems. In the framework of modal superposition method, these fields are the eigenvectors related to the free-vibration problem. If the preliminary problems provide linearly independent fields, then the coefficients of the linear combination are the reduced-state variables of the ROM. On the contrary, one has to build a reduced-basis of the subspace spanned by the known fields. The Proper Orthogonal Decomposition (POD) aimed to create such a basis from a set of time dependent fields.

The POD model reduction method has been used in a wide range of nonlinear incremental simulations and optimization problems in fluid mechanics [START_REF] Daescu | Efficiency of a POD-based reduced second-order adjoint model in 4D-Var data assimilation[END_REF], in materials science [START_REF] Ganapathysubramanian | Design across length scales: a reduced-order model of polycrystal plasticity for the control of microstructure-sensitive material properties[END_REF][START_REF] Monteiro | Computational homogenization for nonlinear conduction in heterogeneous materials using model reduction[END_REF], in thermal science [START_REF] Balima | Comparison between the modal identification method and the POD-Galerkin method for model reduction in nonlinear diffusive systems[END_REF] and real-time surgery simulation [START_REF] Niroormandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF]. This method comes from the Karhunen-Loève expansion [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF][START_REF] Loève | Probability theory[END_REF] developed for statistical analyses. A POD basis is an optimal basis of the state subspace spanned by forecasted states possibly related to different simulations of the response of the detailed model. The optimality of the POD has been stated in [START_REF] Holmes | Low-dimensional models of coherentstructures in turbulence[END_REF]. The snapshot POD proposed by L. Sirovich [START_REF] Sirovich | Turbulence and the dynamics of coherent structures partI : coherent structures[END_REF] had an important contribution to the development of POD basis for the reduction of nonlinear incremental problems. Some details about this method are given in Section 4. During the same period, P. Ladevèze proposed a competitive decomposition method dedicated to elastoplastic models named the radial loading decomposition (RLD) [START_REF] Ladevèze | Sur une famille d'algorithmes en mècanique des structures[END_REF]. This decomposition is performed by the LATIN method [START_REF] Ladevèze | Sur une famille d'algorithmes en mècanique des structures[END_REF] without solving the full incremental detailed model. The LATIN method aims at building a linear problem to forecast the response of the nonlinear detailed model. This linear problem is defined over the entire time interval. Therefore, a time and space separated representation enables to define the decomposition of the forecasted state evolution. Various LATIN algorithms were proposed to solve optimization problems [START_REF] Allix | A new multi-solution approach suitable for structural identification problems, Comput[END_REF][START_REF] Boucard | A multiscale strategy for structural optimization[END_REF] and multiscale problems [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF]. The ***similarities*** between the RLD and the POD have been stated in [START_REF] Ryckelynck | Réduction a priori de modèles thermomécaniques[END_REF].

In practice, during the solution of a series of problems, the POD ROM tends to be less convenient as the current problem differs from the preliminary problems. The main drawback of the classical POD and the snapshot POD methods is the lack of an adaptive procedure in order to modify the ROM basis according to an error indicator. A first adaptive procedure was proposed in [START_REF] Ryckelynck | Réduction a priori de modèles thermomécaniques[END_REF] using the LATIN method. But the non-incremental scheme used in the LATIN to build linear problems does not facilitate the extension of the method to any mechanical constitutive law. An a posteriori error estimation and a basis adaptation have been proposed in [START_REF] Nguyen | A posteriori error estimation and basis adaptivity for reduced-basis approximation of nonaffine-parametrized linear elliptic partial differential equations[END_REF] for the reduction of nonaffine-parametrized linear elliptic partial differential equations. This is an interesting approach since a series of nonlinear problems involves a series of parametrized linear problems. But it is no easy task to extend this method to mechanical problems involving internal variables. In most cases, the time is not a simple parameter of the partial differential equations. These difficulties have led us to adopt an incremental approach. Different incremental adaptive strategies were proposed to adapt the reduced-basis approximation: for nonlinear thermal problem [START_REF] Ryckelynck | A priori hypereduction method: an adaptive approach[END_REF],

for solving some transfer equations [START_REF] Verdon | An adaptive ROM approach for solving transfer equations[END_REF], for kinetic theory models [START_REF] Ammar | On the reduction of kinetic theory models related to finitely extensible dumbbells[END_REF], and for any known state evolution [START_REF] Ryckelynck | On the "a priori" model reduction: Overview and recent developments[END_REF]. These papers, excepted the last one, are related to the APHR method using Krylov subspaces in order to extend the subspace spanned by the ROM basis. The recent advance on Hyper Reduction methods proposed in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables International Journal for Numerical[END_REF] makes possible the extension of the APHR method to the nonlinear mechanical models involving internal variables. But, the expansion using Krylov subspace was not retained. The Krylov approach increases the number of balance residual evaluations. Therefore the computational time devoted to residuals becomes too expensive in case of complex constitutive laws. We propose to extend the ROM basis by using the solution of the FE equations over few time increments. As proposed in [START_REF] Markovinovic | Accelerating iterative solution methods using reduced-order models as solution predictors[END_REF], a ROM predictor accelerates the Newton-Raphson iterations related to the FE equations.

The multi-level APHR method provides state estimations by summing Finite Element corrections to ROM predictions. The reduced-state variables of the ROM being global and the FE variables being nodal unknowns, the proposed approximation is multi-level. If at the end of a time increment, the ROM prediction is accurate enough, no Finite Element correction is performed. On the contrary, the results of the FE correction make possible to adapt the ROM. When a FE correction has been performed an approximate state evolution is known. Therefore we can apply the adaptive algorithm proposed in [START_REF] Ryckelynck | On the "a priori" model reduction: Overview and recent developments[END_REF] dedicated to known state evolutions. The mechanical state is taken into account to expand the subspace spanned by the reduced-basis related to the ROM. To master the growth of the ROM, the adaptive procedure involves a POD of the reduced-state variables. As a result of the Hyper Reduction method [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables International Journal for Numerical[END_REF], a specific spatial integration scheme is introduced when computing the reduced-state variables related to the ROM prediction. The constitutive equations are integrated over a Reduced Integration Domain (RID), provided that this equations are local. The quality of the ROM prediction is evaluated thanks to the norm of a truncated residue of the FE balance equations. The computation of the FE correction can be performed using any classical incremental algorithm. It could be a classical Newton-Raphson algorithm or a Newton-Raphson algorithm coupled with a parallel solver such as the FETI method [START_REF] Farhat | Implicit pa-rallel processing in structural mechanics[END_REF] or a mixed domain decomposition method [START_REF] Ladevèze | On mixed and multiscale domain decompositon method[END_REF].

The paper is organized in the following manner. Section 2 introduces the formulation of the mechanical problem involving internal variables. The formulation of the equations related to the multi-level APHR method is detailed in Section 3. Section 4 reports the numerical results elucidating the usefulness of the proposed method for simplified sensitivity analysis. A modification of a parameter related to the material behaviour and a modification of a geometrical parameter are considered. The results obtained with a sequential solver and those obtained with the FETI parallel solver are discussed in this section. The computational time savings are compared to those obtained with the snapshot POD. Section 5 is the conclusion of this paper. Some indications about the ongoing work are given.

Formulation of the continuous model

The continuous model of concern is a parametrized mechanical model. We denote {p} the column of the model parameters. These parameters can be related to design choices, material coefficients or loading conditions. We consider a series of mechanical problems related to a series of parameter values ({p} α ) α=1...N . ***The various models of the series are supposed to be similar. Therefore, we assume that the potential geometrical modifications can be represented using an unchanged tessellation of the domain with modified positions of vertices.*** We assume that displacements and strains are small. The stress tensor σ is a nonlinear function of the strain history depending on the parameters {p} α :

σ = Σ(ε τ , τ ≤ t; {p} α ) (1)
where Σ is a formal operator that must be defined by constitutive equations.

The continuous medium is occupying a domain Ω({p} α ). The nonlinear system is analyzed over a time interval ]0, T ]. The displacement field at time t is defined on Ω({p} α ) and it is denoted by u(., t, {p} α ).

The boundary ∂Ω({p} α ) of Ω({p} α ) is denoted by ∂ U Ω({p} α ) ∪ ∂ F Ω({p} α ). On ∂ U Ω({p} α ), there is the Dirichlet condition u(., t, {p} α ) = 0 for all t. On ∂ F Ω({p} α ), there is a given force field F(., t, {p} α )

depending on time t and parameters {p} α . The displacement field belongs to a function space V defined by:

V = u(., t, {p} α ) ∈ H 1 (Ω({p} α )) | u | ∂ U Ω({p} α ) = 0 (2)
The statement of the mechanical problem is the following. We want to find an estimation of the displacement field u ∈ V defined by the constitutive equations and the principle of virtual work:

Ω({p} α ) ε(u * ) : Σ(ε(u) , τ ≤ t, {p} α ) dΩ - ∂ F Ω({p} α ) u * . F(x, t, {p} α ) dΓ = 0 ∀ u * ∈ V (3)
where u * is a test function.

In cases of stable materials, the knowledge of the operator Σ is sufficient to find a solution of the mechanical balance equations (3). According to the framework of the irreversible thermodynamic processes, a constitutive law can be defined by a choice of: internal variables z, a free energy w(ε, z, {p} α ) and a pseudo potential of dissipation ϕ * [START_REF] Germain | Continuum Thermodynamics[END_REF]. Some conjugated variables Z are associated to the internal variables z using the definition of the dissipation. The conjugated variables and the internal variables are connected by the following equation of state:

Z = - ∂w ∂z (4)
Complementary constitutive equations can be proposed without introducing a pseudo potential of dissipation, provided that the Clausius Duhem Inequality is fulfilled ( the rate of entropy production must not be negative):

ż = B(Z, {p} α ) (5)
In case of standard formulation of the constitutive equations, the complementary constitutive equations are deduced from the pseudo potential of dissipation ϕ * (Z, {p} α ) such that:

B(Z) = ∂ϕ * ∂Z (6)
The initial state of the material is defined by a given initial condition:

z |t=0 = z ini (7)
The stress σ being one of the conjugated variables, the set of equations ( 4) to [START_REF] Halpen | Sur les materiaux standards generalises[END_REF] defines the operator Σ.

3 The multi-level method

Definition of function spaces

Three functional subspaces of V are introduced to facilitate the formulation of the proposed method. The displacement field is forecasted by a ROM prediction and a FE correction. Due to the proposed adaptive procedure the ROM is modified during the incremental treatment of the series of mechanical problems.

Therefore we introduce the superscript (n) which denotes the version of the ROM. Using the FE method [START_REF] Zienkiewicz | Finite Element Method[END_REF] we obtain the detailed model to be reduced. It involves m h degrees of freedom (q j (t, {p} α )) j=1...m h .

The component q j is a nodal displacement connected to the displacement field by the shape function N j .

The FE model is deduced from the continuous model by using a subspace of V denoted V h such that :

V h = {u ∈ V |∃ {q} ∈ ℜ m h , u(x) = j=m h j=1 N j (x, {p} α ) q j ∀x ∈ Ω({p} α )} (8)
We assume that the topology of the mesh does not depend on the geometrical parameters possibly involved in {p} α . ***To study simple design modifications, only the position of some nodes of the mesh is supposed to change. The list of elements and the list of nodes connected to each element remain unchanged.***

The ROM is defined by a subspace of

V h denoted V (n)
ROM such that:

V (n) ROM = {u ∈ V h |∃ {a} (n) ∈ ℜ s , u(x) = k=s k=1 ψ (n) k (x) a (n) k ∀x ∈ Ω({p} α )} (9)
where (ψ

(n) k ) k=1.
..s are the shape functions of the ROM. They are deduced from the FE shape functions by using a reduction matrix [A]

(n) such that :

ψ (n) k = j=m h j=1 N j (x, {p} α ) A (n) jk ∀x ∈ Ω (10)
As proposed in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables International Journal for Numerical[END_REF] a RID is introduced to define the balance condition related to the ROM according to the Hyper Reduction method. The basic idea is to select only few nodal balance equations of the FE model in order to state the formulation of the balance condition of the reduced-order model. Therefore only few elements contribute to this formulation. These elements are forming the RID. The most similar idea found in the literature concerns the reduction of chemistry models [START_REF] Banerjee | An adaptive reduction scheme to model reactive flow[END_REF], despite it is not applied on FE models. A mathematical programming based approach provides a selection a few species and chemical reactions in order to reduce the model.

Let's denote [Π] the rectangular matrix selecting the FE balance equations such that Π ij is equal to one if the degree of freedom number j is the i th selected degree of freedom. Else, Π ij is equal to zero. A truncated test function u * Π can be associated to any test function u * = N. {q} * such that

u * Π = N. [Π] T . [Π] . {q} * (11)
Therefore the following property is fulfilled:

if u * ∈ V h then u * Π = N. [Π] T . [Π] . {q} * ∈ V h . The RID Ω Π ({p} α
) is defined as the support of all the truncated test fields:

Ω({p} α ) = Ω Π ({p} α ) ⊕ Ω Π ({p} α ) (12
)

Ω Π ({p} α ) = x ∈ Ω({p} α ) | ∀ {q} * ∈ ℜ n , u * Π (x) + ε(u * Π )(x) = 0 (13)
where

||ε(u * Π )|| 2 = ε(u * Π ) : ε(u * Π ).
The hyper reduced-balance condition is obtained using the space function V

(n)

Π such that:

V (n) Π = { u ∈ V (n) ROM |∃ {a} (n) ∈ ℜ s , u(x) = j=n j=1 N j (x) q j ∀x ∈ Ω, with {q} = [Π] T . [Π] . [A] (n) . {a} (n) } (14)
The hyper reduced-balance condition is the following:

ΩΠ({p} α ) ε(u * ) : Σ(ε(u) , τ ≤ t, {p} α ) dΩ - ∂ F ΩΠ({p} α ) u * . F(x, t, {p} α ) dΓ = 0 ∀ u * ∈ V (n) Π (15)
where

∂ F Ω Π ({p} α ) = ∂ F Ω({p} α ) ∩ Ω Π ({p} α )
. This equation involves s scalar equations, where s is the size of the ROM. To obtain a well-posed reduced-problem, the number of selected FE balance equations must be sufficient. It can't be lower than the size of the reduced-basis. A methodology to build this selection of equations is proposed in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables International Journal for Numerical[END_REF]. Obviously, if all the degrees of freedom are selected to build the RID then the RID covers all the domain Ω({p} α ). It is clear that the local computations related to such balance condition are restricted to the RID. The equation [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF] enables to forecast the reduced-state variables related to the displacement field. The ROM shape functions being global, the displacement field is globally forecasted. But to forecast the entire mechanical state, we must extend the estimated internal variables from the RID to the entire domain. As proposed in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables International Journal for Numerical[END_REF] the APHR ROM have two bases, one related to the displacements and one related to the internal variables. The second basis is denoted

(Υ (n) k ) k=1...ξ .
The column of the reduced-internal variables is denoted {b} (n) . The internal variables are extrapolated from the RID to the full domain by fitting the reduced-internal variables inside the RID. The detailed equations of this extrapolation procedure are given in the next section.

The multi-level incremental algorithm

The multi-level approximation is the following:

u(x, t, {p} α ) = u (n) ROM (x, t, {p} α ) + δu h (x, t, {p} α ) u (n) ROM ∈ V (n) ROM δu h ∈ V h (16) z(x, t, {p} α ) = z ROM (x, t, {p} α ) + δz(x, t, {p} α ) (17)
with

z ROM (x, t, {p} α ) = k=ξ k=1 Υ (n) k b (n) k (t, {p} α ) (18)
Using a numerical one-step time integration scheme, one can forecast different states of the system at different instants. According to an incremental formulation, the mechanical state is assumed to be known at time t i . The unknowns are the state variables at time t i+1 . Several stages are introduced to forecast the mechanical state over the time increment ]t i , t i+1 ]:

• stage one : the reduced-state variables related to the displacement field are forecast with δu h = 0 and δz = 0

• stage two : the internal variables are extended using the basis related to these variables

• stage three : an error indicator related to the accuracy of the ROM is evaluated

• stage four : if the ROM prediction is not accurate enough, then the FE correction (δu h , δz) is computed

• stage five : if the FE correction has been performed, then the reduced-bases are adapted by using the results of the correction stage (stage four).

If the ROM prediction is accurate enough, then no FE computation is performed. In practice, the adaptation of the ROM is achieved only using few FE increments. Moreover, in case of nonlinear equations or parallel computing, the iterative solver takes advantage of the prediction provided by the ROM. The error indicator is the norm of a truncated residue of the FE balance equations. This is the residue of the FE equations selected for the creation of the RID. At the end of stage four, a full state estimation is known. Then the adaptive algorithm [START_REF] Ryckelynck | On the "a priori" model reduction: Overview and recent developments[END_REF] devoted to known state evolutions can be applied.

The adaptation stage involves an expansion of the subspaces related to the ROM and a selection of the most significant events using a POD of the reduced-state variables. The extension of the bases is provided by orthogonal contributions extracted from the state corrections. During the adaptation of the ROM the reduced-state variables related to the previous computations are updated. A norm of the load magnitude is introduced to define a relative error indicator such that:

F 2 = ∂ F ΩΠ({p} α ) F(x, t, {p} α ) . F(x, t, {p} α ) dΓ (19)
The detailed formulation of the five stages of the algorithm is given below.

Stage one, the unknown are the reduced-state variables {a} (n) (t i+1 , {p} α ) and the internal variables z ROM (x, t i+1 , {p} α ) over Ω Π ({p} α ) such that:

ΩΠ({p} α ) ε(u * ) : Σ(ε(u (n) ROM ) , τ ≤ t, {p} α ) dΩ - ∂ F ΩΠ({p} α ) u * . F(x, t, {p} α ) dΓ = 0 ∀ u * ∈ V (n) Π (20)
Stage two, the unknown are the reduced-internal variables {b} (n) (t i+1 , {p} α ) and the internal variables

z ROM (x, t i+1 , {p} α ) localized in Ω Π ({p} α ) such that: {b} (n) (t i+1 , {p} α ) = arg min {y} H({y}) (21)
H({y}) = ΩΠ({p} α ) z ROM (x, t i+1 , {p} α ) - k=ξ k=1 Υ (n) k (x) y k 2 z dΩ (22)
z ROM (x, t i+1 , {p} α ) = k=ξ k=1 Υ (n) k (x) b (n) k (t i+1 , {p} α ) ∀x ∈ Ω Π ({p} α ) (23)
Stage three, evaluation of the error indicator η ROM related to the truncated residue {R}:

R j = ΩΠ({p} α ) ε(u * ) : Σ(ε(u (n) ROM ) , τ ≤ t, {p} α ) dΩ - ∂ F ΩΠ({p} α ) u * . F(x, t, {p} α ) dΓ (24)
with

u * = N j (x, {p} α ) (max ζ Π ζj ) (25)
η ROM = {R} (26)
Stage four, if η ROM < ǫ R F then δu h (x, t i+1 , {p} α ) = 0 and δz(x, t i+1 , {p} α ) = 0, else these corrections are such that:

Ω({p} α ) ε(u * ) : Σ(ε(u (n) ROM + δu h ) , τ ≤ t, {p} α ) dΩ - ∂ F Ω({p} α ) u * . F(x, t, {p} α ) dΓ = 0 ∀ u * ∈ V h (27) Stage five, if η ROM ≥ ǫ R then the subspace V (n)
ROM can be adapted by using δu h . We only consider the residue related to the orthogonal projection of δu h into V

ROM . This residue is denoted δ ⊥ u h and it is such that:

δ ⊥ u h (x, t i+1 , {p} α ) = δu h (x, t i+1 , {p} α ) -δu h (x, t i+1 , {p} α ) ( 28
)
δu h (x, t i+1 , {p} α ) = k=s k=1 ψ (n) k (x) δa (n) k (t i+1 , {p} α ) (29)
{δa} (n) (t i+1 , {p} α ) = arg min {y} Ω({p} α ) δu h (x, t i+1 , {p} α ) - k=s k=1 ψ (n) k (x) y k 2 dΩ (30)
Then, an extended basis ψ

(n+1/2) k k=1...s+1
is built preserving the previous predictions such that:

ψ (n+1/2) k = ψ (n) k k ≤ s (31) ψ (n+1/2) s+1 = 1 δ ⊥ u h δ ⊥ u h (32) a (n+1/2) k (τ, {p} α ) = a (n) k (τ, {p} α ) τ ≤ t i+1 k ≤ s (33)
a (n+1/2) k (τ, {p} β ) = a (n) k (τ, {p} β ) τ ≤ T k ≤ s β < α (34)
a (n+1/2) s+1 (τ, {p} β ) = 0 τ ≤ T β < α (35)
a

(n+1/2) s+1 (τ, {p} α ) = 0 τ ≤ t i (36)
a (n+1/2) s+1 (t i+1 , {p} α ) = δ ⊥ u h (37)
A POD decomposition of the reduced-variables is performed to avoid a constant growth of the size of the ROM. This POD consists in finding the vectors ({V } l ) l=1...s+1 maximizing the following projection on the forecast reduced-state variables:

λ (n+1) l = β=α-1 β=1 T 0 {a} (n+1/2)T (t, {p} β ) . {V } l 2 dt + tj+1 0 {a} (n+1/2)T (t, {p} α ) . {V } l 2 dt {V } l 2 (38)
within the conditions {V } l = 1 and λ l ≥ λ l+1 . The vectors ({V } l ) l=1...s+1 are the eigenvectors of the covariance matrix [C] (n+1/2) such that:

[C] (n+1/2) = β=α-1 β=1 T 0 {a} (n+1/2) (t, {p} β ) . {a} (n+1/2)T (t, {p} β ) dt + tj+1 0 {a} (n+1/2) (t, {p} α ) . {a} (n+1/2)T (t, {p} α ) dt (39)
The covariance matrix takes into account the previous results to preserve the ability of the ROM to model the related mechanical states. Then, the main events are selected using the following criteria:

[V ] = {V } 1 , ..., {V } s with λ s > ǫ P OD λ 1 and λ s+1 ≤ ǫ P OD λ 1 (40)
s is the new size of the ROM. At last, the new ROM is a POD reduction of the previous ROM:

ψ (n+1) l = k=s+1 k=1 ψ (n+1/2) k V kl ∀l = 1... s (41) a (n+1) l (t, {p} β ) = k=s+1 k=1 a (n+1/2) k (t, {p} β ) V kl ∀l = 1... s (42)
A similar adaptive procedure is performed to build (Υ

(n+1) k , b (n+1) k
) k=1...ξ using δz h (x, t i+1 , {p} α ). This adaptation ends the treatment of the time increment ]t i , t i+1 ].

At the end of each simulation, a last POD of the reduced-state variables is performed. The parameters of the method are the RID and the coefficients ǫ R and ǫ P OD . The initial ROM ((ψ

(0) k , a (0) k ) k=1...s , (Υ (0) k , b (0)
k ) k=1...ξ ) can be empty (s = 0 and ξ = 0). In such a case, the first increment of the simulation starts with one full FE increment. The adaptive algorithm proposed in [START_REF] Ryckelynck | On the "a priori" model reduction: Overview and recent developments[END_REF] being used, the following property is fullfiled. In case of full integration (Ω Π = Ω), the decomposition obtained at the end of the incremental computation is the POD of a state prediction having an accuracy related to ǫ R . The proof of this property can be found in [START_REF] Ryckelynck | On the "a priori" model reduction: Overview and recent developments[END_REF].

Parallel computing

The correction stage is a classic FE computation using an initial guess state. Therefore a classical domain decomposition method can be used to solve this problem on a parallel computer. We chose to apply the FETI method [START_REF] Farhat | Implicit pa-rallel processing in structural mechanics[END_REF]. But any mixed method is also convenient [START_REF] Ladevèze | On mixed and multiscale domain decompositon method[END_REF]. The solution of stage four is provided by a FETI solver coupled to a Newton-Raphson algorithm [START_REF] Feyel | Computation of a Triaxial Specimen Using a Polycrystalline Model[END_REF][START_REF] Feyel | Application du calcul parallèle aux modèles à grand nombre de variables internes[END_REF].

A full parallel version of the multi-level APHR method is achieved by using the same domain decomposition for all the stages. Each domain of the decomposition has its own contribution to the RID and its own part of the global shape functions of the ROM. We use the single-program multiple-data approach. Therefore each processor is forecasting the global reduced-state variables. But the computations of contributions to reduced-residuals and to reduced-tangent stiffness matrices are distributed over all the processors. Each processor is summing the contributions of all the other processors. The smaller the ROM size the smaller is the amount of exchange between processors during stage one. Moreover, the POD of the reduced-state variables is duplicated on each processor without any exchange of data between them.

***This method was implemented in the existing Finite Element Analysis software Z-set (www.nwnumerics.com). We modified the linear solver and we extended the data structure of the

The elastoviscoplastic constitutive law

This section aimed to compare the multi-level APHR method, the snapshot POD method, and the FE method. A classic and simple elastoviscoplastic example is considered in order to facilitate the reproduction of the presented results by the reader. The strain tensor is split into the elastic strain tensor and the plastic strain tensor:

ε = ε p + ε e (43)
The accumulated plastic strain p is defined by:

p = t 0 2 3 εp τ : εp τ dτ (44)
Two internal variables homogenous in strain are used: the plastic strain tensor ε p and a tensor ξ related to a kinematic hardening such that:

z = ε p ξ , z 2 z = ε p : ε p + ξ : ξ (45)
The free energy is assumed to be of the following form (isotropic elasticity):

w(ε, ε p , ξ) = 1 2 E 1 + ν (ε -ε p ) : (ε -ε p ) + ν E (1 -2ν)(1 + ν) ((ε -ε p ) : I) 2 + 1 2 c ξ : ξ (46)
in which E is the Young modulus, ν is the Poisson ratio and c is a scalar coefficient. c is the first parameter of the mechanical model. The conjugated variables associated to the internal variables are the stress tensor and the kinematic hardening tensor Ξ such that:

σ = - ∂w ∂ε p , Ξ = - ∂w ∂ξ (47)
We then obtain the following state laws:

σ = E 1 + ν ε e + ν E (1 -2ν)(1 + ν) (ε e : I) I , Ξ = -c ξ (48)
The yield surface is defined by the function f :

f (σ, Ξ) = (σ D -Ξ) : (σ D -Ξ) 1/2 -G o (49)
in which σ D is the deviatoric part of the stress tensor: σ D = σ -1 3 (ε e : I) I. The complementary constitutive equations are provided by a standard formulation. The pseudo potential of dissipation is such that:

ϕ * (σ, Ξ) = K θ + 1 f (σ, Ξ) K θ+1 + (50)
Hence, we obtain the normality rule:

εp = ∂ϕ * ∂σ = ṗ σ D -Ξ ((σ D -Ξ) : (σ D -Ξ)) 1/2 (51) ξ = ∂ϕ * ∂Ξ = -εp (52) ṗ = f (σ, Ξ) K θ + ( 53
)
The initial conditions are:

ξ |t=0 = ε p |t=0 = 0.

The series of mechanical problems

The system of interest for sequential computations, is a square plate (80 mm x 80 mm x 5 mm, E = 98000 A series of three simulations is considered. The variable parameters of the models are the radius r and the material coefficient c:

MPa, ν = 0.3, G o = 330
{p} T = {r, c}
The three simulations are related to {p}

d i < r 1 d i + ∆d i = 1 r + 1 r i - 1 r 1 r2 -1 r 1 r1 -1 r (54) if d i ≥ r ∆d i = 0 (55)
where d i is the distance between the closest hole center and the node number i and r = 40mm ***( r must be greater than the hole radius and smaller than the distance between the hole and the closest external boundary of the domain). Using the mapping defined by equations ( 54) and (55), the smaller the radius modification the smaller the distances between the new Gauss points and the old Gauss points are. The meshes have always the same number of Gauss points and the same number of nodes. Therefore the transfer of the reduced-bases is straightforward. The reduction matrices related to the reduced-basis remain unchanged. We just take into account node position modifications.***

The maximum accumulated plastic strain is reached at point Q shown on Figure 1 for PM. The σ 22ε 22 curve related to the point Q is shown on Figure 1 for PM. An implicit Euler scheme and a Newton-Raphson algorithm have been chosen for the solution of the nonlinear time dependent equations.

The time interval is regularly split into 151 time steps. The SM model involves 3,200 linear hexahedral elements, 332,800 scalar internal variables and 19,926 degrees of freedom related to nodal displacements.

The PM model involves 51,200 linear hexahedral elements, 5,324,800 scalar internal variables and 312,942 degrees of freedom related to nodal displacements. Four subdomains have been used to perform the domain decomposition related to the FETI method. These domain decomposition is shown on Figure 1.

The simulations have been performed using the ZeBuLon code [START_REF] Besson | Large scale object-oriented finite element code design[END_REF] and the very efficient DSCPACK linear solver [START_REF] Raghavan | Domain-Separator Codes for the parallel solution of sparse linear systems[END_REF] for sequential simulations.

Usual strategy using the snapshot POD

The first FE simulation related to {p} 1 is used has preliminary problem in order to create the snapshots.

These snapshots are the displacements forecast by the FE model at the end of each time step. According to the snapshot POD method the POD basis (φ k ) k=1...γ is such that:

φ k (x) = i u(x, t i , {p} 1 ) d ik (56)
{d} k maximize µ k (57)

µ k = T 0 Ω u(x, t, {p} 1 ).φ k dΩ 2 dt Ω φ k . φ k dΩ (58)
The key point about the snapshot POD is the reduction of the size of the eigenproblem related to the maximization of µ k by using a linear combination of the snapshots. The POD basis is used to performed the reduced-simulations related to the parameters {p} 2 and {p} 3 . A fixed functional space is therefore introduced to define the POD ROM. This subspace is denoted V P OD . It is such that:

V P OD = {u ∈ V h |∃ {a} ∈ ℜ s , u(x) = k=s k=1 φ k (x) a k ∀x ∈ Ω({p} α )} (59)
The reduced-state variables {a} and approximated FE internal variables are defined by the constitutive equations and a classical Galerkin procedure:

Ω({p} α ) ε(u * ) : Σ(ε(u) , τ ≤ t, {p} α ) dΩ - ∂ F Ω({p} α ) u * . F(x, t, {p} α ) dΓ = 0 (60) ∀ u * ∈ V P OD with u ∈ V P OD (61) for α ∈ {2, 3} (62)
The constitutive equation being solved over the entire domain, no reduced-basis is required to forecast the internal variables. A prediction of these variables is provided at each integration point of the mesh.

Comparison of five strategies

The purpose of this section is the evaluation of the capability of the ROMs to forecast the modification of The size of the ROMs are reported in Table 1. These results are very similar. During the third simulation (α = 3), the ROM adaptations produced by APS and AAPHRS increase the size of the ROM. The ROM version correspond to the accumulated amount of FE increments used to extend the APHR ROM. The amount of FE increments related to each series of simulations is shown in Table 2. Each FE increment involves several solution of linear FE problem. The amount of full size linear problems which have been solved is reported in Table 3. Thanks to the initial prevision provided by stage one of the multi-level APHR method the savings related to the amount of global linear problems is bigger than the savings related to the amount of full FE increments. One can observe in Table 2 that only few FE increments were used to adapt the APHR ROMs.

The five strategies are compared considering various computational times and the accuracy of the maximal accumulated plastic strain. Each series of simulations is performed using the same computer in case of sequential computation or the same set of processors in case of parallel computations. The duration of parallel simulations is given by the processor having the longest computational time. The quality of each simulation is reported in Table 4. The computational time related to each simulation and the total computational time related to each strategy are reported in Table 5. The computational time savings for PM are in good agreement with savings related to the amount of global linear problems. It is not the case for SM because local computations have a greater impact on computational time in case of small size problem.

The Hyper Reduction method has a bigger contribution to the computational time savings in such a case.

Let's consider the full series of simulations. Computational time saving must be considered depending on the expected accuracy of the simulations. If the most accurate results are expected then the FE strategy must be used. But APS is also very convenient in such a case and less expansive than FES. Due to the adaptive procedure, APS is the most accurate strategy using a ROM. If the required level of error is 30%, then all the results can be considered. In such a case, the fastest strategy is FAPHRS. But if the required level of error is lower than 10%, then PODS and FAPHRS are failing to provide convenient results for the third simulation (α = 3) of the series. In such a case the computational savings provided by PODS and FAPHRS can't be considered. Therefore, the fastest strategy is AAPHRS.

The Hyper Reduction method produces error amplifications. As mentioned in [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables International Journal for Numerical[END_REF], the bigger the RID the smaller are the error amplifications. By increasing the size of the RID the AAPHRS strategy can tends to APS strategy. By decreasing the parameter ǫ R the FAPHR strategy can tends to APHRS.

Let's consider now each simulation one by one. The maximum computational time saving is 76 % for SM and 97 % for PM. These results are provided by the multi-level APHR methods in case of FAPHRS.

It turns out that the best results concerning the accuracy or concerning the computational time are provided by the multi-level APHR method. APS is less efficient than AAPHRS and FAPHRS. The Hyper reduction technique is very useful to perform fast simulations. The proposed method enables us to choose between accurate simulations or fast simulations. Moreover, one can observe that during the first simulation (α = 1) APS, AAPHRS and FAPHRS are more efficient than the FES while the PODS is necessary less efficient than FES. This makes the multi-level APHR method very efficient for parallel computing.

Conclusion

No preliminary finite element solution is needed to build APHR bases. The multi-level APHR method can be considered as a new efficient solver. Even without knowing an initial ROM, the multi-level APHR method is more efficient than FE method. The extension of this solver to parallel or sequential iterative solver is straightforward. This extension has been performed using a parallel FETI solver. Obviously the better initial ROM the faster is the multi-level APHR simulation. Accurate simulations can be performed without using Hyper Reduction. The Hyper Reduction of the equations provides significant computational time savings. Fast simulations can be performed using Hyper Reduction without ROM adaptation. The multi-level APHR method enables to choose gradually between fast or accurate simulations by changing two inputs: the RID and the parameter ǫ R . The best computational saving shown in this paper is 97%. This method is very efficient in case of series of simulations related to modifications of the mechanical model.

***We addressed in this paper simple geometrical modifications by moving some nodes of the mesh. More complex geometrical changes involving remeshing techniques can be certainly taken into account provided that mesh transfer operations are applied to the reduced-basis. The better the reduced-bases are transferred from the previous mesh to the new mesh, the better the quality of the initial ROM will be.

  nonlinear solver ZeBuLon. The new data structure involves the reduction matrix [A] (n) the reduced state variables {a} (n) and the index of the ROM version n for the reduced-bases related to displacements and internal variables. The reduction matrix contains nodal values for displacement basis vectors. It contains Gauss point values for internal variable basis vectors. The reduced integration scheme is activated by changing the list of the elements of the mesh.*** 4 Numerical results

  MPa) with a hole in the middle (radius r). The system of interest for parallel computations, is a square plate (320 mm x 320 mm x 5 mm, E = 98000 MPa, ν = 0.3, G o = 330 MPa) having 16 holes in it (radius r). In the text below, we refer to the first model as sequential model (SM) and to the second one as the parallel model (PM). Due to the symmetry of both domains and boundary conditions, only one 8 th of the plates are modelled. The boundary conditions are shown on Figure1.

T 1 = 2 = 3 =

 123 {r 1 , c 1 }, {p} T {r 1 , c 2 }, {p} T {r 2 , c 1 } using r 1 = 20 mm, r 2 = 22 mm, c 1 = 35000 MPa, c 2 = 28000 MPa. ***To represent the modification of the holes radius by preserving the elements of the mesh, the radial position of the nodes of the mesh is changed such that:*** if

  the maximum accumulated plastic strain p max at point Q. The reference values are provided by the FE simulations for each α in {1, 2, 3}. Therefore the reference strategy consists in three FE simulations. The second strategy is related to the snapshot POD. It involves a preliminary FE simulation and two reduced-simulations using a fixed POD basis as presented in the previous section. The third strategy introduces the adaptive procedure without Hyper Reduction. The three simulations (α ∈ {1, 2, 3}) are performed using the multi-level APHR method with Ω Π = Ω, ǫ P OD = 10 -7 , ǫ R = 0.01 for the first simulation (α = 1) and ǫ R = 0.1 for the last simulations (α ∈ {2, 3}). The fourth strategy introduces the Hyper Reduction method. The third strategy is restarted by changing the reduced-integration domain Ω Π . This domain is fixed during the treatment of the three simulations (α ∈ {1, 2, 3}). The RIDs related to the SM and PM models are shown on Figure1. Ω Π involves 1,083 and 16,640 elements respectively for SM and PM. The fifth strategy is reproducing the fourth using less accurate criteria to adapt the ROM. The parameter ǫ R have been set equal to 1 during the last simulations (α ∈ {2, 3}). These five strategies are denoted respectively FE strategy (FES), POD strategy (PODS), A Priori strategy (APS), accurate APHR strategy (AAPHRS) and fast APHR strategy (FAPHRS).

6 Acknowledgments

 6 Anyway, as mentioned above, if the quality of the transferred ROM (the initial ROM of the new simulation) is too bad, the proposed adaptive algorithm will add the convenient corrections during the new simulation.*** The ongoing works are related to the reduction of models involving damage growth, finite strains and life duration predictions. The case of large number of simulations involved in the series of problems needs deeper research. The extent of the efficiency of the multi-level APHR method remains an open question.***The APHR method should be very efficient to simulate smooth nonlinear state evolutions. For instance, in case of contact problems (with or without friction), if the contact zones jump from a place to an other at each time increment and if this process is very sensitive to the parameters of the model (stochastic evolution), we can't expect the APHR method to be efficient. But if the contact areas have a smooth evolution in time, the APHR method should be efficient.*** This work is related to the MELOXEL project partially funded by the French National Research Agency (ANR).

Figure 1 :

 1 Figure 1: SM and PM models

Table 1 :

 1 Size of ROMs

		α = 1 α = 2 α = 3
	for SM
	PODS		6	6
	APS	6	6	7
	AAPHRS	6	6	8
	FAPHRS	6	6	6
	for PM
	PODS		7	7
	APS	6	6	9
	AAPHRS	6	7	7
	FAPHRS	6	6	6
		20

Table 2 :

 2 Amount of full FE increments

		α = 1 α = 2 α = 3 Total Saving
	for SM
	FES	151	151	151	453	0%
	PODS	151	0	0	151	66%
	APS	78	7	33	118	74%
	AAPHRS	67	0	24	91	80%
	FAPHRS	67	0	0	67	85%
	for PM
	FES	151	151	151	453	0%
	PODS	151	0	0	0	66%
	APS	78	9	29	116	74%
	AAPHRS	54	0	3	57	87%
	FAPHRS	54	0	0	54	88%

Table 3 :

 3 Amount of global linear problems

		α = 1 α = 2 α = 3 Total Saving
	for SM
	FES	256	257	292	805	0%
	PODS	256	0	0	256	68%
	APS	78	9	58	145	82%
	AAPHRS	71	0	55	126	84%
	FAPHRS	71	0	0	71	91%
	for PM
	FES	248	247	282	777	0%
	PODS	248	0	0	248	68%
	APS	78	11	51	140	82%
	AAPHRS	63	0	4	67	91%
	FAPHRS	63	0	0	63	92%

Table 4 :

 4 Relative error between the approximate and the FE maximal accumulated plastic strain

		α = 1	α = 2	α = 3	maximum
	for SM
	PODS < 0.5% < 0.5 %	26%	26%
	APS < 0.5% < 0.5 % < 0.5 %	< 0.5 %
	AAPHRS < 0.5% < 0.5 %	9 %	9 %
	FAPHRS < 0.5% < 0.5 %	29 %	29 %
	for PM
	PODS	2%	< 1 %	17%	17%
	APS	2%	< 1 %	< 1 %	2 %
	AAPHRS	2%	3 %	6 %	6 %
	FAPHRS	2%	4 %	26 %	26 %
			21

Table 5 :

 5 Computational times

		α = 1	α = 2	α = 3	Total	Total CPU Saving Maximum CPU Saving
	for SM
	FES	601 s	603 s	719 s	1923 s	0%	0%
	PODS	671 s	362 s	335 s	1368 s	29%	53%
	APS	583 s	389 s	589 s	1561 s	19%	35%
	AAPHRS	419 s	143 s	476 s	1038 s	46%	76%
	FAPHRS	418 s	144 s	168 s	1038 s	62%	76%
	for PM
	FES 27 067 s 25 740 s 28 661 s 81 468 s	0%	0%
	PODS 27 587 s	1 524 s	1 542 s	30 655 s	62%	95%
	APS	9 914 s	2 685 s	6 878 s	19 477 s	76%	90%
	AAPHRS 10 376 s	751 s	2 995 s	14 112 s	83%	97%
	FAPHRS 10 376 s	724 s	732 s	11 832 s	85%	97%