B Bensaid

X Boddaert

P Benaben

R Gwoziecki

R Coppard

Reliability of OTFTs on flexible substrate: Mechanical stress effect

Flexibility will significantly expand the application scope of electronics, particularly large-area electronics. Over the last ten years, printed organic electronic is believed to be one of the next major technological breakthroughs in the field of microelectronic and the use of printing technology to process organic fieldeffect transistors (OFETs) opens promising perspectives for low cost, large area circuits integrated on flexible, plastic substrates. With amorphous polymer-based thin films transistors acceptable electrical performances are now achieved with relatively good stability at ambient air. In literature a lot of work has been devoted to study degradation of device characteristics under bias stress conditions but only few papers deal with the mechanical behaviour. In this paper, we review our first reliability results obtained on flexible organic thin film transistors under mechanical stresses. The variations of electrical characteristics under bending tests, both in compression and tension, have been studied. Using specific equipment, we have also evaluated the reliability of transistors under cyclic bending tests. The stress dependency of the transfer characteristic deviates from the one observed for inorganic material like silicon.

Introduction

Organic transistors based on conjugated polymers or small molecules have been intensively studied in recent decades, and their performance has continually improved [1] with understanding and progress on chemistry design and synthesis of solution-processed organic semiconductors [2] [3]. Organic thin films transistors (OTFTs) have the advantages of flexibility, large area capability, low cost, and low temperature processing, and so have various potential applications such as in displays [4], radio frequency identification tags [5], and large area sensors [START_REF] Someya | Proc. Natl. Acad. Sci[END_REF] on flexible substrates. These applications are commonly subjected to bending strains, which could modify the performance of organic devices. One of the next challenges is now to demonstrate the robustness and reliability of OTFTs to support the use in the frame of these applications. While degradation of device characteristics under bias stress conditions has been extensively studied [START_REF] Sirringhaus | [END_REF], only few papers deal with the effect of mechanical stress on the electrical characteristics. In particular, small molecules-based OTFT and changes in mobility observed in bent organic TFT devices have been reported to result from changes in the energy barrier to hole hopping due to variation of the intermolecular spacing under a bending strain [8]. In this paper, we aim to investigate the effect of mechanical stresses on our printed OTFT characteristics.

Material and methods

In this work, the mechanical behaviour of a printed p-channel organic field effect transistor is studied. Samples are achieved in the CEA-Liten laboratories. The structure of the transistor is a "top gate -bottom contacts" design (Fig. 1.b). The semi-conducting layer is a screen-printed triarylamine copolymer film, which is a well known family of air stable and solution-processed P-channel semiconductor [9]. The nonplanar conformation of the triarylamine and the lack of order make this class of semi-conducting material suitable for printing and annealing processes. Gold source and drain contacts were patterned by laser ablation on a 125 µm plastic PolyEthylene Naphthalate (PEN) substrate from Dupont Teijin sputtered with 30 nm of gold. Above the organic semiconductor, a 1µm thick low k perfluoropolymer dielectric layer (Cytop ® ) [10] and a silver conductive gate have been also screen printed. Each 5 cm x 5 cm test sheet contains 30 transistors with channel lengths of 50 µm, 80 µm and 100 µm and channel widths ranging from 500 µm to 2 mm. For fatigue tests, larger test sheets 8 cm x 20 cm, also composed 30 transistors have been designed to fit the machine. Electrical characterizations of the OTFTs have been performed in the dark, in air without encapsulation, using a Keithley 4200 semiconductor analyser combined with a semi-automatic Süss Microtech prober. Transistors parameters such as saturation mobility, threshold voltage and I on/off ratio were extracted following the 1620TM IEEE standard test procedure. Figure 1.a shows the typical output characteristics without stress, source-drain current I DS is monitored as a function of source-drain voltage V DS and gate voltage V GS was changed from 0 to -40 V in -10 V steps. The measured mobility µ was 0.03 cm 2 /V.s and the on/off ratio was 10 5 .

Results and discussion

Bending tests both in tension and compression have been performed while recording the electrical characteristics of the transistor (Fig. 2). Curvature radii R from 50 mm down to 3 mm have been tested in tension. Due to probing problems, only large radii between 50 and 30 mm have been achievable in compression mode. Since the transistor layers are very thin (µm or below) compared to the PEN substrate thickness t s , absolute value of the mechanical strain 1 in the transistor can be evaluated from bending radius using the following formula (Eq.1) [11]. By convention, compressive strains are negative values while tensile strains are positive.

|1| = (t s / 2R) (1)

Transistors are tested with two different stress directions under identical conditions so that their source-drain current paths may be arranged precisely perpendicular and parallel to the direction of strain, referred to as perpendicular position and parallel position, respectively. The shape of I DS -V GS and I DS -V DS characteristics under bending strain remains close to the initial one (Fig. 3.a and Fig. 3.b respectively). The transistor still works normally under bending stress but presents electrical performance variations. Gold source and drain electrodes give good ohmic contact for all I DS -V DS curves and contact resistance does not change drastically. Leakage current increases by one decade : from 1pA to 10 pA. The electrical properties of transistors were also measured while changing bending radius. Normalized saturation current and mobility have been plotted versus the mechanical strain 1 with the parallel position (Fig. 4.a and Fig. 4.b respectively). Saturation current increases clearly with increasing tensile strain. The change in saturation current is +12% at bending radius R=3 mm. On the opposite, in compression mode, although we have only investigated small strains, the saturation current tends to decrease. Such large changes cannot be explained only by the deformation of the device structures at R=3 mm, corresponding to strain of only 2.08%. Indeed, using the classical Poisson model (Eq.2) [START_REF] Kittel | Solid State Physics[END_REF] where 1 P and 1 i are Poisson coefficients of semiconductor and gate insulator, respectively, changes in the structural parameters of OTFT, namely, L (channel length), W (channel width), and t i (the thickness of gate insulator) are evaluated and their insertion in the equation Eq.3 shows that for 2.08% of strain geometrical deformation causes only about -2% of saturation current.

W 1 W (1 -2 P . 1) ; L 1 L (1+ 1) ; t i 1 t i (1 -2 i . 1) (2)
Finally when the stress is released, we recover the initial electrical characteristics of the transistor. The change in source-drain current is reversible and reproducible. This reversibility shows that this phenomenon is not associated with the time-dependent drift of electronic performance.

We evaluated field-effect mobility µ using the equation (Eq.3) [START_REF] Kimissis | Organic Field Effect Transistors[END_REF], where I DS is the source drain current, C i is the insulator capacitance, and V T is the threshold voltage, in the saturation mode and plotted it normalized by the value without strain µ 0 in the Fig. 3. In this figure the mobility increases by 8% when strain is 2.08% of tension. The observed strain effect on the mobility of these organic TFTs with amorphous polymer semiconductor is the opposite of the results obtained with crystalline or polycrystalline organic semiconductors like pentacene films [START_REF] Yang | [END_REF] [15] when mobility decreases with increasing tensile strain and increases with increasing compressive strain. In this last case, because of the size and the organization of pentacene small molecules, mobility change by mechanical stress is quite reasonable since the transport properties of organic molecular systems follow the hopping model rather than the band model. Under tensile strain, as described in a theoretical study [16] energy barrier for hole hopping increases due to the larger spacing between ordered molecules, resulting in an decreasing in the mobility. In contrast, in the case of our study, the polytriarylamine has an amorphous and non-planar structure with a very long polymer chains muddled randomly [17] which can explain the opposed behaviour compared to pentacene molecules. Interestingly, we note that parallel position and perpendicular position show similar strain dependences of mobility, although the directions of strain are perpendicular and parallel to the direction of the source-drain current path for these two different positions. Such an isotropic electrical property indicates that transport is also consistent with the hopping model in amorphous thin films, in which a conducting path from the source electrode to the drain electrode is formed by coupled grain chains oriented in random directions [START_REF] Kimissis | Organic Field Effect Transistors[END_REF].

I D =(µC i /2) (W/L) (V G -V T )² (3)
We have also studied the total resistance variation versus the mechanical strain. Using the transmission line method (TLM), we have observed that the channel resistance decreases while the source and drain contact resistance increases under a tensile bending strain. At tensile strain of about 2% which is close to the maximum strain of gold film [18], the device is still functional and the contact resistance variation shows that charge injection still occurs at source/semiconductor interface; confirmed by the I DS -V DS characteristic. This result could be explained by the fact that, in the transistor configuration, gold electrodes are encapsulated by the organic semiconductor layer which minimizes the degradation of gold/semiconductor interface. The channel resistance change is consistent with the mobility variation under tensile stress. Thus, the increase of the mobility could be explained by a modification of the polymer structure. Indeed tensile stress could induce a slip mechanism of the polymer chains which improve the relative order in the polymer material [19] which leads to an increase of the conductivity [20]. Electrical observation and parameters variation seem to confirm the coherence of this explanation based on literature analysis and although we have no enough clear proof yet of this structural modifications under strain. So the next step of this study will be a detailed structural analysis under stress by AFM measurements and RAMAN spectroscopy [21] [22] in order to highlight the structural variation of the polytriarylamine during mechanical stress. Complementary to these static tests and in order to investigate the recovery performance after stressing transistors, some preliminary fatigue tests are performed using a dedicated cyclic bending machine (Fig. 5).

A lot of bending cycles with bending radius of 5 mm, which corresponds to a strain of 1.25%, are applied to the OTFTs. First, the transistors sheet is bent from flat R=1 to R=5 mm, and immediately released up to flat state at the highest frequency of our system (15 cycles per minute). Measurements are performed at 50, 100, 300, 500, and 1000 bending cycles. The normalized saturation current as a function of the number bending cycles is shown in Fig. 6. We observed an unexpected increase of the saturation current I Dsat by ~10% as soon as we reached 50 cycles. I Dsat continues increasing with the number of cycles. After 1000 times, the augmentation is estimated to 60%. 4 weeks after this experiment, we measured again our devices and observed a change of the electrical characteristic; indeed, saturation current has decreased from 60% to 20% (Fig. 5). That means that after stress, a slow relaxation phenomenon occurs. In spite of the lack of catastrophic degradation after a lot of cycles, we have no good stability under mechanical fatigue contrary to encapsulated organic transistors which are more stable under a lot of bending cycles [23]. But this performance remains better than amorphous-silicon field effect transistors (a-Si FETs) which fail mechanically by periodic cracks in the amorphous layers of the TFT at a strain of about 0.5% [24] [25]. More experiments are necessary with increasing the number of cycles and investigating the relaxation after fatigue to get a better understanding of this phenomenon which could be a real problem for the functionality of electronic systems.

Conclusion

In summary, we have studied the robustness of organic field-effect transistors (OFET's) made of amorphous triarylamine copolymer semiconductor on a PEN flexible substrate under mechanical stress. Under tensile bending tests, we have observed a reversible increase of the saturation current and mobility of about 10% for a mechanical tensile strain of 2%. This could be explained by an improvement of the relative order of the semiconductor structure. Further ways of analysis could still be investigated in order to achieve more information about the structural behaviour of the amorphous active layers. Indeed, using micro-Raman spectroscopy which is a sensitive technique of structure study, polymer organisation generated by mechanical strain can be brought out after analysing Raman spectra. Besides that, the AFM microscopy, with phase image in tapping mode, can supply the illustration of these structure variations due to mechanical strain. We have also performed preliminary cyclic bending tests and observed a 60% increase of the saturation current after 1000 cycles with a slow relaxation phenomenon. Different transistor geometries and orientations will be tested to get a better understanding of these mechanisms.

Fig. 1 :

 1 Fig.1: (a) illustration of "top gate -bottom contacts" organic TFT on PEN substrate. (b) Typical dc characteristics without stress, showing source-drain current I D as function of source-drain voltage V D . Gate voltage V G is changed from 0 to -40 V in -10 V steps.

Fig. 2 :

 2 Fig.2: Photograph and Schematic illustration of electrical characterization by tips probing under tensile bending. Using halfcylinder forms with Curvature radiuses R from 50 mm (1 = 0.13%) down to 3 mm (1 = 2.08%).

Fig. 3 :

 3 Fig.3: I DS -V GS characteristics at V DS = -40V (a) and I DS -V DS characteristics at V GS = -40V (b) with bending strains 1 from 0 (flat position) to +2.08% (bending radius of 3 mm).

Fig. 4 :

 4 Fig.4: Relative change in (a) saturation current I Dsat and (b) mobility µ plotted as function of strain tensile strain 1. Compressive strains are negative values (only 1 from -0.13% to -0.21%) while tensile strains are positive (1 from +0.13% to +2.08%).

Fig. 5 :Fig. 6 :

 56 Fig.5: bending fatigue machine. Different curvature radii can be used (5 mm, 10 mm and 25 mm) with a maximum speed of 15 cycles by min.

Acknowledgments

This work is performed in the frame of the collaboration between CEA-LTEN of Grenoble and ENSM of St-Etienne. The authors want to thank these two institutions and their heads for supporting and financing this research activity.


David Uystepruyst
email: david.uystepruyst@univ-valenciennes.fr

Mame William-Louis

Emmanuel Creusé

Serge Nicaise

François Monnoyer

Efficient 3D numerical prediction of the pressure wave generated by high-speed trains entering tunnels

Keywords: Computational fluid dynamics, Euler equations, Three-dimensional simulation, Cartesian mesh, Sliding grid, Non-reflecting boundary conditions, High-speed trains

The topic of this paper is to present a new methodology for the threedimensional numerical simulation of the entrance of high-speed trains in a tunnel. The movement of the train is made thanks to a technique of sliding meshes and a conservative treatment of the faces between two domains. All parts of the development are thought with the aim to reduce the computational time. In particular, non reflecting boundary conditions for non-structured three-dimensional meshes are developed in order to limit the calculation domain. Validations of the methodology are presented on different test cases.

Introduction

During the train entrance in a tunnel, the air initially inside the tunnel is strongly compressed by a piston effect in the train head neighborhood and gives rise to a compression wave that propagates at the sound speed in the tunnel. During its propagation in the tunnel, this wave can be damped in the case when the railway is on ballasted track, Ozawa [START_REF] Ozawa | Effect of ballast on pressure wave propagating through tunnel[END_REF], or can be steepened in the case when the railway is on concrete slab track, due to the nonlinear effects [START_REF] Matschke | Full scale tests on pressure wave effects in tunnels[END_REF]. On reaching the distant tunnel exit, the compression wave is partially reflected back into the tunnel as a rarefaction wave, and a part of it emerges as a pressure pulse, called micro-wave. In certain circumstances this micro-wave may be strong enough to generate annoying sonic disturbances. Indeed, it can generate a booming noise up to 140-150 db or more, Maeda [START_REF] Maeda | Effect of shape of train nose on compression wave generated by train entering tunnel[END_REF]. The magnitude and the duration of this wave are strongly connected with the temporal gradient of the initial compression wave. Hence it seems necessary to simulate with accuracy the generation of the compression wave both in amplitude and in temporal gradient. It is a challenging task, even for the current most advanced computers systems.

The flow generated by a train which circulates in a tunnel is threedimensional, non stationary, turbulent and compressible. A complete study has to take into account all these aspects. However, comparing the results obtained by Ogawa and Fujii [START_REF] Ogawa | Numerical investigation of three-dimensional compressible flows induced by a train moving into a tunnel[END_REF] or by Shin and Park [START_REF] Shin | Numerical study of flow characteristics of the high speed train entering into a tunnel[END_REF] that used the Navier-Stokes equations with the works of Réty and Grégoire [START_REF] Réty | Numerical simulation of the pressure wave generated by a train enters a tunnel[END_REF][START_REF] Réty | Numerical investigation of tunnels extensions attenuating the pressure gradient generated by a train entering a tunnel[END_REF], Mok and Yoo [START_REF] Mok | Numerical study on high speed train and tunnel hood interaction[END_REF] as well as Yoon [START_REF] Yoon | Prediction and validation on the sonic boom by a high-speed train entering a tunnel[END_REF], that have chosen an Eulerian approach, one realizes that the initial wave is perfectly simulated by both approaches. Actually, the pressure jump caused by the nose train entry does not depend on the viscosity [START_REF] Hara | Aerodynamics problem when a train running into a tunnel with a large velocity[END_REF], contrary to the second jump for which the wall friction effects are involved. However, this second jump, quasi-linear, does not modify the maximum pressure gradient. That is why Euler equations can be used to model this phenomenon.

During the displacement of the train, a relative movement with respect to the tunnel and to the ground can be distinguished. For the management of this movement, several techniques can be used:

• The chimera method, used by Ogawa [START_REF] Ogawa | Numerical investigation of three-dimensional compressible flows induced by a train moving into a tunnel[END_REF] or Yoon [START_REF] Yoon | Prediction and validation on the sonic boom by a high-speed train entering a tunnel[END_REF] in a simplified version, also called overlapping meshes method. The computations are independently made in several subdomains and an interpolation has to be done for the transfer of the variables between the domains.

To limit the interpolation errors, it is necessary to refine the common zones and consequently the computational cost can become very heavy, especially in the case of explicit time-integration schemes for which stability criteria induce small time steps.

• The method of meshes deformation, used by Réty [START_REF] Réty | Numerical simulation of the pressure wave generated by a train enters a tunnel[END_REF][START_REF] Réty | Numerical investigation of tunnels extensions attenuating the pressure gradient generated by a train entering a tunnel[END_REF] or Mok [START_REF] Mok | Numerical study on high speed train and tunnel hood interaction[END_REF] consists in building a new and more regular mesh during the simulation when the elements of the mesh become too stretched due to the deformations. An interpolation is also necessary to make the link between the old mesh and the new one.

• Finally, the technique of sliding meshes consists, as for the chimera method, to subdivide the computational domain into several subdomains, but without any intersection. This method is well suitable for a rectilinear movement like the one studied in this paper. Shin and Park [START_REF] Shin | Numerical study of flow characteristics of the high speed train entering into a tunnel[END_REF] use this technique with an interpolation for the communication between domains.

During the simulations, the flow around the train is supposed to be stationary before it enters the tunnel. To accelerate the convergence towards this stationary state, non reflecting boundary conditions are needed to allow the reduction of the computational domain and to eliminate numerical disturbances.

This paper presents the development of a three-dimensional simulation based on the resolution of the Euler equations. The first originality of this work is to consider the displacement of the train by using a sliding meshes technique with a conservative treatment of the faces between two domains.

The second one is to extend the non reflecting boundary conditions of Giles [START_REF] Giles | Non-reflecting boundary conditions for Euler equation calculations[END_REF] to three-dimensional and unstructured meshes in this sliding meshes framework for general artificial boundaries. These contributions allow to get more accurate as well as low expensive numerical simulations, by reducing the computational domain without spurious non-physical reflection.

The paper is organized as follows: section 2 describes the methodology of sliding meshes. Then, section 3 is devoted to the model, the numerical scheme used, the non-reflecting boundary conditions implementation and their validation on an overpressure/depression test case. Finally, some validations are performed on real cases to evaluate the jump of the pressure and of its temporal gradient induced by the entrance of a train in a tunnel. Here, experimental results [START_REF] Maeda | Effect of shape of train nose on compression wave generated by train entering tunnel[END_REF][START_REF] Réty | Numerical simulation of the pressure wave generated by a train enters a tunnel[END_REF] confirm the reliability of the numerical methodology used in the present study. The methodology [START_REF] Waymel | Numerical simulation of the unsteady three-dimensional flow in confined domains crossed by moving bodies[END_REF] adopted here consists in subdividing the domain into two (non-overlapping) sub-domains. The first one (sliding domain part) contains the train and can be set in motion, as shown on Figure 1. The second domain encapsules the first one. It contains walls of the tunnel and the external domain. The sliding domain contains two additional parts: the "upstream zone" and the "downstream zone" whose limits are fixed. As the elements of the mesh arrive at the end of the "upstream zone", due to the movement of the domain, they are automatically re-injected at the beginning of the "downstream zone". With this approach, there is no need to regenerate the mesh and the length of the sliding domain can be kept constant. Both sub-domains are independently meshed by means of an automatic Cartesian mesh generator [START_REF] Deister | Self-organizing hybrid cartesian grid generation and application to external and internal flows, Numerical Flow Simulation III[END_REF] based on a boundary triangular mesh. The grid generator relies directly on an octree structure. The boundary cells are truncated and a merging procedure allows to avoid the convergence problem caused by small cells (see [START_REF] Deister | Self-organizing hybrid cartesian grid generation and application to external and internal flows, Numerical Flow Simulation III[END_REF] for all details). Both sub-domains are connected to each other at common interfaces where a reconstruction of the cells faces is made on each side to preserve the conservative properties of the numerical scheme.

Methodology of sliding meshes

Equations and numerical scheme

Equations

As mentioned in the introduction, the pressure wave jump generated by the train entry in the tunnel is modeled by the 3D Euler equations, corresponding to the conservation of mass, quantity of movement and energy :

∂ t U + ∇ • H(U) = 0, (1)
where U is the vector of conservative variables and H(U) the flux tensor:

U =            ρ ρu ρv ρw ρE            , H(U) =            ρu ρv ρw ρu 2 + p ρuv ρuw ρvu ρv 2 + p ρvw ρwu ρwv ρw 2 + p ρuh 0 ρvh 0 + pv ρwh 0            .
Here ρ stands for the volumic mass, u, v, w are the three components of the speed vector V = (u, v, w) t and E is the energy per unit of mass.

Under the perfect gas assumption, the static pressure p is as usual given by:

p = (γ -1)ρ E - 1 2 (u 2 + v 2 + w 2 ) ,
where γ = 1.4 and h 0 is the total enthalpy defined by:

h 0 = E + p ρ .
Now, let us introduce the translation velocity of the sliding mesh V t = (u t , v t , w t ) t , and the relative velocity of the fluid V r = (u r , v r , w r ) t = V-V t .

In the sliding domain, equation (1) becomes:

∂ t U + ∇ • H (r) (U) = 0. (2)
Where the fluxes tensor H (r) (U) reads:

H (r) (U) =            ρu r ρv r ρw r ρuu r + p ρuv r ρuw r ρvu r ρvv r + p ρvw r ρwu r ρwv r ρww r + p ρu r h 0 + pu t ρv r h 0 + pv t ρw r h 0 + pw t            .

Numerical scheme

Integrating equation (1) on a control volume V, we get

∂ t V U dV + S H n (U) dS = 0,
where the volumic integral of the divergence of H becomes a boundary integral due to Green-Ostrogradski formula, where H n (U) = H(U) • n, n being the normal vector of the control volume.

After discretization of the computational domain, Roe's scheme [START_REF] Roe | Approximate riemann solvers, parameter vectors, and difference schemes[END_REF] is used to determine the numerical approximation of the flux H n (U) at the interfaces of the control volumes. This numerical flux can be defined as the sum of a centered contribution and of a correction which depends on the eigenvectors and the eigenvalues of the Jacobian matrix ∂H n (U)/∂U. These last ones are well known, see Hirsch [START_REF] Hirsch | Numerical Computation of Internal and external Flows[END_REF], in the case of the classical Euler equations.

In order to calculate the numerical flux in the sliding domain, we have to re-determine them in the case of the Euler equations with a relative velocity.

Let us denote D n the Jacobian matrix ∂H n (U)/∂U and D n (r) the Jacobian matrix ∂H n (r) (U)/∂U. Clearly:

H n (r) (U) = H n (U) -(V t ) n U , so that: D n (r) = D n -(V t ) n I ,
where I is the identity matrix and

(V t ) n = V t • n.
Consequently, the eigenvalues, λ r) are obtained by subtracting (V t ) n to the eigenvalues of D n :

(j) r (1 ≤ j ≤ 5), of D n (
λ (1) r = λ (2) r = λ (3) r = (V r ) n , λ (4) r = (V r ) n + c and λ (5) r = (V r ) n -c, (3)
where c = γp/ρ is the sound speed.

The eigenvectors, useful for the construction of Roe's scheme, remain the same.

With this formulation that takes into account the relative velocity, we can apply second order Roe's scheme for the spatial discretization. The temporal integration is realized with a predictor-corrector scheme of Van Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme v.a. second order sequel to godunov's method[END_REF]. To preserve the TVD property of the scheme during the passage at the second order, the limiter for non-structured mesh of type Barth and Jespersen [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] is used. In order to determine the fluxes between the two meshes, a calculation of the common interface (Figure 2), is made first, using the Sutherland-Hodgman [START_REF] Sutherland | Reentrant polygon clipping[END_REF] algorithm for convex polygon clipping. Figure 3 The translation velocity of the sliding mesh can a priori be considered as an additional difficulty. Nevertheless, by definition of the sliding meshes technique, the normal vectors of the interfaces are always orthogonal to the translation vector of the mesh, and thus

V t • n = 0.
So the translation velocity does not interfer in the calculation of the flux at the interfaces. Hence a classical flux computation can be performed there. This allows a conservative calculation of the fluxes at the interfaces between 10 the two meshes. The validation of the method is performed on the shock tube defined on Two configurations have been considered, as shown on Figure 5. The first one consists in a single domain, while the second one is subdivided into two domains. In this configuration, the second domain can be set in motion.

In order to measure the capability of the method to handle the propagation of strong gradients across the sliding interface, a numerical sensor is located in the fixed domain just behind the sliding domain, represented by a cross into Cartesian meshes on Figure 6. On this figure, the sliding domain is slightly longer than the fixed domain to ensure the cyclic regeneration of cells. Qualitatively, it can be shown that the conservative method is in very good accordance with the reference solution obtained with the single block calculation, both v equal 0 and v equal 70 m/s. In order to quantify these results the L ∞ norm of non-dimensional pressure errors are summarized in Table 1.

When the sliding domain does not move, the obtained solution has to be the same than the solution of the single block. The existing small error is due to the fact that volumic meshes are slightly different between both It can be showed that the sliding grid method gives the same convergence order than the single block. Moreover, the sliding domain curve gets closer to the single block curve while the space mesh size decreases.

p * -(p * ) ref ∞ p * -(p * ) ref ∞ / (p * ) ref ∞ v=0 m/s 9.

Non reflecting boundary conditions

The non reflecting boundary conditions allow to make the calculations without the appearance of non physical reflections. Giles [START_REF] Giles | Non-reflecting boundary conditions for Euler equation calculations[END_REF] introduced such a type of boundary conditions for the two-dimensional Euler equations and in the case of boundary faces parallel to one axis of the domain. These conditions are based on a Fourier decomposition of the solution at the boundary. Then Nicoud extended Giles' conditions for arbitrary boundary faces in the case of a subsonic 2D [START_REF] Nicoud | Defining wave amplitude in characteristic boundary conditions[END_REF] and 3D [START_REF] Nicoud | On the amplitude of waves in characteristic boundary conditions for Euler conditions[END_REF] outflow.

The boundary condition presented here is defined for an arbitrary threedimensional face and applied to outflow as well as inflow boundaries. Furthermore, the translation velocity of the sliding domain is taken into account.

Practically, let us suppose that we have to define the fluxes through a face centered in F (see Figure 9) belonging to a boundary on which some non reflecting boundary conditions are imposed (see for instance the boundary of the external domain in dotted line on Figure 1). For that purpose, Euler equations are written in a characteristic formulation. First, let us define the three fluxes F (r) (U), G (r) (U) and K (r) (U), respectively, as the first, the second and the third column of the flux tensor H (r) (U). Equation (2) can be rewritten:

∂ t U + A (r) ∂ x U + B (r) ∂ y U + C (r) ∂ z U = 0, (4)
where A (r) , B (r) , C (r) are the corresponding Jacobian matrices diagonalized, with the eigenvalues given by (3):

∂F (r) (U) ∂U , ∂G (
Λ (r) = LD n (r) L -1 = diag(λ (1) r , λ (2) r , λ (3) r , λ (4) r , λ (5) r ).
Characteristic variables are defined by the relation valid for arbitrary variations δ, see Hirsch [START_REF] Hirsch | Numerical Computation of Internal and external Flows[END_REF]:

δW = LδU, δU = L -1 δW.
Premultiplying (4) by L, we obtain the five equations:

                     ∂ t W 1 + V r .∇W 1 = 0 ∂ t W 2 + V r .∇W 2 + 1 2 cs 1 .(∇W 4 + ∇W 5 ) = 0 ∂ t W 3 + V r .∇W 3 + 1 2 cs 2 .(∇W 4 + ∇W 5 ) = 0 ∂ t W 4 + (V r + cn).∇W 4 + c(s 1 .∇W 2 + s 2 .∇W 3 ) = 0 ∂ t W 5 + (V r -cn).∇W 5 + c(s 1 .∇W 2 + s 2 .∇W 3 ) = 0,
where (n, s 1 , s 2 ) forms an orthonormal basis, see Figure 9. In a full expansion, it is equivalent to

                     ∂ t W 1 + (V r ) n ∂ n W 1 + (V r ) s 1 ∂ s 1 W 1 + (V r ) s 2 ∂ s 2 W 1 = 0 ∂ t W 2 + (V r ) n ∂ n W 2 + (V r ) s 1 ∂ s 1 W 2 + (V r ) s 2 ∂ s 2 W 2 + 1 2 c(∂ s 1 W 4 + ∂ s 1 W 5 ) = 0 ∂ t W 3 + (V r ) n ∂ n W 3 + (V r ) s 1 ∂ s 1 W 3 + (V r ) s 2 ∂ s 2 W 3 + 1 2 c(∂ s 2 W 4 + ∂ s 2 W 5 ) = 0 ∂ t W 4 + ((V r ) n + c)∂ n W 4 + (V r ) s 1 ∂ s 1 W 4 + (V r ) s 2 ∂ s 2 W 4 + c(∂ s 1 W 2 + ∂ s 2 W 3 ) = 0 ∂ t W 5 + ((V r ) n -c)∂ n W 5 + (V r ) s 1 ∂ s 1 W 5 + (V r ) s 2 ∂ s 2 W 5 + c(∂ s 1 W 2 + ∂ s 2 W 3 ) = 0. (5)
The outflow condition of Giles, extended by Nicoud [START_REF] Nicoud | On the amplitude of waves in characteristic boundary conditions for Euler conditions[END_REF] for a 3D arbitrary face is:

∂ t W 5 = -(V r ) s 1 ∂ s 1 W 5 -(V r ) s 2 ∂ s 2 W 5 -(V r ) n (∂ s 1 W 2 + ∂ s 2 W 3 ) . ( 6
)
Using the fifth equation of ( 5), condition ( 6) can be rewritten in its spatial form:

∂ n W 5 = ∂ s 1 W 2 + ∂ s 2 W 3 . (7)
With the strategy used by Nicoud, we extend the original Giles' conditions to an arbitrary three dimensional boundary face. Results obtained are given in Table 2. Those conditions are now named NRBC (Non Reflecting Boundary Conditions). Those expressions allow us to make the extrapolation discussed below.

In order to define the fluxes at the boundary, U F is defined from U P using a suitable extrapolation (see Figure 9), that reads:

U F = U P + (∇U) P • PF = U P + (∂ n U) P PF • n + (∂ s 1 U) P PF • s 1 + (∂ s 2 U) P PF • s 2 . Flow Conditions outgoing ∂ n W 5 = ∂ s 1 W 2 + ∂ s 2 W 3 ∂ n W 1 = 0 incoming ∂ n W 2 = ∂ s 1 W 2 ∂ n W 3 = ∂ s 2 W 3 ∂ n W 4 = - 1 2 (∂ s 1 W 2 + ∂ s 2 W 3 )
= L -1 (∂ n W) P .
Two cases have to be considered:

• Outflow: Only the fifth characteristic variable leaves the domain: ∂ n W j

(1 ≤ j ≤ 4) are evaluated by a discretization inside the computational domain and ∂ n W 5 is derived from ∂ s 1 W 2 + ∂ s 2 W 3 using the Giles conditions (Table 2).

• Inflow: We proceed similarly using this Table 2 to express ∂ n W j (1 ≤ j ≤ 4) and ∂ n W 5 is discretized inside the domain.

To validate these non reflecting boundary conditions, a domain D 0 defined as the unit cube is considered. This domain is discretized by 100 elements in each direction. As initial condition, a slightly de-centered Gaussian function of pressure is imposed:

                             p 0 =          p ∞ 1 + µe -1/( √ x 2 +(y-0.15) 2 +z 2 -0.3) 2 for x 2 + (y -0.15) 2 + z 2 < 0.3, p ∞ else, u 0 = v 0 = w 0 = 0 m.s -1 , T 0 = 293 K, ρ 0 = p 0 R T 0 , (8)
where p ∞ = 101325 P a.

Equation ( 1) is solved with this initial condition.

This initialization condition implies an overpressure followed by a depression. Hence arrived at the limit of the domain, the flow is at first outgoing and then ingoing. This allows to validate both conditions given in Table 2. In order to evaluate the results obtained on this domain D 0 , a reference solution was computed using a larger domain which can be considered as "infinite", since the waves do not reach its boudaries during the simulation time involved.

Graphs from Figure 10 show the propagation of the wave generated by the initial data [START_REF] Mok | Numerical study on high speed train and tunnel hood interaction[END_REF]. The reference value is represented by the black continuous lines, while the value given by the calculation onto the truncated domain is It can be seen that the extrapolation method implies an important drift.

                       ρ F = ρ P + -→ ∇ρ P • -→ P F u F = u P + -→ ∇u P • -→ P F v F = v P + -→ ∇v P • -→ P F w F = w P + -→ ∇w P • -→ P F p F = p P + -→ ∇p P • -→ P F ,
Table 3 shows the quantitative results. The error generated by the NRBC method is about four times lower than the error of the extrapolation method. Finally, as noticed on the various graphs and on the error on Table 4 indicates requisite CPU time for those simulations. The calculation on the domain D 0 , with Giles' method, is slightly more time consuming than the extrapolation method on this same domain. Nevertheless, the calculation on the "infinite" domain, which is two times larger than D 0 in each direction, requires about 8 times more computational efforts. With the degree of convenience between results obtained on domain D 0 with Giles' method and on the "infinite" domain, it is possible to consider domain D 0 rather than the "infinite" domain and, then, to strongly reduce the computational time.

p * -(p * ) ref ∞ / (p * ) ref ∞ NRBC method Extrapolation method x = z = 0 1.92 × 10 -2 7.86 × 10 -2 x = z = 0.25 5.71 × 10 -2 2.31 × 10 -1
In order to realize a convergence analysis, the reference solution obtained on the "infinite" domain with a mesh size of ∆x = 1/100 m is compared with numerical solutions resulting of calculations performed on the domain

D 0 with mesh size ∆x = 1/(20 * i) m for i = 1, . . . , 4.
Results are displayed on Figure 11. This graph gives the relative L ∞ norm error of the non-dimensional pressure time history for the sensor coordinates (0.25;0.5;0.25), located at the boundary. It can be shown the error decreases with the mesh size.

Solid walls

In order to consider the correction of the flux of the numerical scheme for faces situated on the boundary of the domain, a ghost cell is built in which the flow velocity V ghost is determined thanks to the flow speed in boundary cell V b :

(V ghost ) t = (V b ) t and (V ghost ) n = -(V b ) n
The other variables are identical. In that way the impermeability condition is satisfied while introducing the flux correction of Roe.

High speed train validations

Reduced scale tests

To study the influence of the nose shape of the train on the initial compression wave, Maeda [START_REF] Maeda | Effect of shape of train nose on compression wave generated by train entering tunnel[END_REF] made measurements on a reduced scale model. The body of the reduced train is cylindrical and three shapes of the nose are considered: conical, parabolic and elliptic.

As stated on Figure 12 of pressure and of its temporal gradient recorded at the sensor are compared with the experimental results on the graphs of Figure 15.

Graphs as well as maxima presented in Table 5 show a very good concordance between the experimental and numerical results. The errors on the jumps of pressure do not exceed 6%, while the errors on the maxima of gradients of pressure are lower than 4%. We notice that the committed errors are more important for the conical profile than for two other cases. The geometrical discontinuity of the conical nose, between nose and body, implies this most important error.

Although the shape of the nose has little influence on the jump of pressure, it

Full scale test

The full scale test considered here is based on the experiment realized in situ by the SNCF (French National Railway Company) [START_REF] Réty | Numerical simulation of the pressure wave generated by a train enters a tunnel[END_REF]. This experiment The interest of this experiment is that our results can be compared with the experimental data but also with the numerical results obtained by the SNCF using the software Tg-flo (three-dimensional Eulerian code).

Tg-flo numerical setup

The mesh is composed of tetrahedra and is unstructured. An arbitrary Lagrange-Euler formulation, which allows to take into account the relative movement between the tunnel and the train, is used to solve the Euler equations using a second order Runge-Kutta scheme for the temporal discretization, and a finite element method for the spatial one.

All the numerical data are summarized in the Table 6.

Tg

Results

The graphs of Figure 19 The determination of the pressure amplitude is not easy from the experimental data. Indeed, as previously mentioned, the pressure signal is decomposed into two successive jumps: the first is consecutive to the train entry, and the second depends on the wall friction effects. Only the first is of interest. SNCF researchers defined the pressure jump value from the time corresponds to the moment where the gradient becomes constant, i.e. slightly after 0.2 s. Actually, this value overestimates the real value and this overestimation can be several decades of Pascals. The computation of the SNCF code slightly underestimates this value of the maximum of pressure, while our code overestimates it. On the other hand the maximum of gradient is strongly overestimated by Tg-flo (11%) and Metas1 (14%). This more important error, for our calculation, can be explained by the difference of the space mesh size.

Metas2 implies an underestimation of the pressure amplitude, but the error on the gradient of pressure is divided by more than 3 to falls to 4% and becomes clearly better than the result obtained by Tg-flo code. As the graph of the gradient shows, the curve given by the second simulation follows the shape of the curve of the experimental solution up to the end of the simulation.

Conclusion

An Eulerian code on Cartesian meshes was developed to simulate the entrance of a train into a tunnel. The technique of sliding mesh is adopted with a conservative treatment at the faces between the sliding and fixed meshes.

To accelerate the convergence towards a stationary state of the flow before the train goes into the tunnel, non-reflective boundary conditions have been implemented. This allowed to reduce the domain of calculation and also to prevent from the reflection of numerical disturbances.

The method is validated on model tests as well as on measurements in situ.

The results show that this approach is well adapted for the simulation of the entrance of a train into a tunnel. With this methodology validated, parametrical studies are in progress to study the influence of a hood adjunction on the temporal gradient of the initial wave.

Figure 1 :

 1 Figure 1: Top view of the computational domain.

3. 3 .Figure 2 :

 32 Figure 2: Construction of common interfaces.

  depicted an example of two polygons clipping. Starting with one edge of the first one, the separating line passing by its vertices is defined to subdivide the plane of the polygon into the inside half-plane containing the first polygon and the remaining outside half plane. The result polygon becomes the intersection between the second polygon and the inside half-plane. This process is repeated for all edges of the first polygon.

Figure 3 :

 3 Figure 3: Common interface construction.

Figure 4 :

 4 Figure 4: Shock tube definition.

Figure 4 .

 4 Figure 4. It is a 4 m length tube separated into two parts by a diaphragm in its middle, with the closed side initially set at a pressure of 1.25 bar. The other part is opened to the free atmosphere.

Figure 5 :

 5 Figure 5: Single block model (left) and configuration with sliding block (right).

Figure 6 :

 6 Figure 6: Sketch of shock tubes Cartesian grids.

  configurations. At the speed of 70 m/s, the error slightly increases. Since the timestep used for the calculations is the same for both configurations, the CFL number differs due to the translation velocity of the sliding block.However, the relative error remains only about 1%. Now, let us consider a convergence analysis. Without any available analytical solution, a calculation is performed on the single block with a very fine mesh size ∆x = 1/32 m which leads to 6 553 600 cells. Then, seven calculations are made with space mesh size ∆x = 1/2(i+1) m for i = 1, . . . , 7 for the single block and for the domain with sliding block, with the translation velocity equal to 70 m/s. The evolution of the L ∞ norm of the non-dimensional pressure error is depicted on the Figure8versus the space mesh size.

Figure 8 :

 8 Figure 8: Convergence analysis of the sliding mesh validation.

2 Figure 9 :

 29 Figure 9: Boundary cell.

  represented by dotted lines blue and circles. The isolines of non-dimensional pressure p * = (p -p ∞ ) /p ∞ on the plane x = 0 (top) as well as the difference of pressure on the lines x = z = 0 (bottom right) and x = z = 0, 25 m (bottom left) are shown at time t = 0, 0015 s. The Giles' formulation is compared with a simple extrapolation method:

  which consists in defining the variables values at the face center F from the variables values at the cell center P and their gradients, all of them being evaluated by their values inside the computational domain. Results of this method are represented by green lines and triangles. The two lines allow to show the efficiency of the boundary conditions of NRBC in the case of an almost one-dimensional exit, where the wave is almost parallel to the boundary: right x = z = 0 m ; and a multidimensional exit, where the wave arrives with a non zero angle of attack: x = z = 0, 25 m.

x = z = 0 Figure 10 :Figure 11 :

 01011 Figure 10: Comparison of NRBC method and extrapolation method with reference data. (a) Isocontours pressure: ( ) numerical results on "infinite" domain, ( ) numerical results of NRBC method on domain D 0 , ( ) numerical results of extrapolation method on domain D 0 .

  , these three forms imply a different variation of section. The tunnel is a cylindrical tube with one sensor of pressure located at 1 m from the entrance. The blocking ratio σ = S tr /S tun , defined as the maximum cross-sectional area of the train to the one of the tunnel, is 0.116 and the train is launched at 232 km/h, which corresponds for an ambient temperature of 20 • C at a Mach number of 0.188. Due to the symmetry of the configuration, only the quarter of the domain has to be considered. So the number of cells is about 300 000. Figure13represents the triangular surfacic meshes of the three configurations, and Figure14represents the corresponding volumic cartesian meshes. The numerical results of the jump

Figure 12 :

 12 Figure 12: Geometrical ratios and form of the nose.

Figure 14 :

 14 Figure 13: Surfacic meshes.

Figure 16 :Figure 17 :

 1617 Figure 16: Comparison between the numerical results of the three nose shapes.

  consists in the entrance of the French Hight Speed Train (TGV) into the tunnel of Villejuste situated on the high-speed southwest line, in France. The train measures 74 m and runs at a speed of 220 km.h -1 before entering in a tunnel of section of area of 46 m 2 , that gives a blocking ratio of σ = S tr /S tun = 0.21.

Figure 18 :

 18 Figure 18: Volumic and surfacic meshes of the TGV.

Figure 19 :

 19 Figure 19: History of the pressure and temporal pressure gradients.

Figure 20 displays

 20 Figure 20 displays pressure isocontours on the train and on the tunnel wall. When the train is at the tunnel entry, the air inside the tunnel is at rest again. After a distance of 15 m, the generation and the propagation of the compression wave are clearly highlighted. The little rectangular piece in front of the tunnel is the "upstream zone", previously mentioned, needful for the cyclic regeneration of cells.

Figure 20 :

 20 Figure 20: Generation and propagation of the compression wave (in Pa).

Table 1 :

 1 Infinite norms of the non-dimensional pressure errors, sliding mesh validation.

		87 × 10 -6	1.29 × 10 -3
	v=70 m/s	7.83 × 10 -5	1.03 × 10 -2

Table 2 :

 2 Spatial formulations of NRBC conditions.

(∂ s 1 U) P and (∂ s 2 U) P are evaluated by a discretization of the gradient inside the domain, (∂ n U) P is evaluated by the knowledge of the spatial variation of the characteristic variables from (∂ n U) P

Table 3 :

 3 Relative error norms of non-dimensional pressure, NRBC validation.

Table 3 ,

 3 values obtained on the truncated domain with the NRBC method are in good accordance with the reference values and no significant reflection can

	be seen.
		domain D 0	domain D 0	infinite domain
		extrapolation method NRBC method
	CPU time (s)	1152	1195	9433

Table 4 :

 4 Calculations CPU time (in second).

Table 5 :

 5 Pressure and its gradient maxima and relative error.

		800				300
	∆P (pa)	200 400 600				100 200	pressure gradient (kpa/s)
		0.005 0	0.01	0.005	0.01	0
		t (s)		t (s)
				(a) circular cone
		800				300
	∆P (pa)	200 400 600				100 200	pressure gradient (kpa/s)
		0.005 0	0.01	0.005	0.01	0
		t (s)		t (s)
			(b) paraboloid of revolution
		800				300
	∆P (pa)	200 400 600				100 200	pressure gradient (kpa/s)
		0.005 0	0.01	0.005	0.01	0
		t (s)		t (s)
			(c) ellipsoid of revolution
	Figure 15: Comparison of jumps (left) and gradients (right) of pressure for the 3 shapes
	of the nose; • • •: experimental results, : numerical results.

Table 6 :

 6 Data relative to the numerical experiments, full scale test.

		-flo (SNCF)	Metas1	Metas2
	domain	half-domain	full domain full domain
	dist. train-tunnel at t=0	≃15 m	160 m	160 m
	nbr of cells	200 000	500 000	1 360 000
	(∆x) min	0.1 m	0.2 m	0.1 m

  compare the different results for the jump of pressure and the temporal gradient of pressure, obtained by the experiment result, the numerical solution of the SNCF code and finally Metas simulations. The values of the maxima are presented in Table7.

		∆P max		∂P/∂t max
		value (pa) diff./ref. (%) value (pa/s) diff./ref. (%)
	experimental	1296		8400
	Tgflo	1291	-1	9330	+11
	Metas1	1332	+3	9600	+14
	Metas2	1251	-3	8795	+4

Table 7 :

 7 Pressure and its gradient maxima and relative error.

is noticed that the paraboloid profile clearly minimizes the maximum of the pressure gradient (see Figure 16). This is due to the fact that the variation of the cross-sectional area of the parabolic nose is constant (see the variations of the areas of the sections of shape of the three configurations represented on Figure 17(b)).

The compression wave is delayed, it appears by chronological order: ellipse, parabola, cone, as shown on Figure 16. For a given value of x on the nose, the section of shape of maximal area is the one of the elliptic nose, as shown on Figure 17. At this point, the pressure is thus more important for the ellipse. This value of pressure is reached later for the parabola, and even later for the cone. This explain the temporal shift of the ascent for the pressure gradient for the three configurations.

A complete simulation is subdivided in two phases. The first one consists in the steady state calculation, with a motionless train. The second one corresponds to the train entry in the tunnel. At the beginning of the displacement, the train is located 15 m downstream of the tunnel.

In order to make this second phase, a mesh deformation is used and a remeshing is done when elements become too small or too stretched.

For the transient simulation (second phase), the computational domain is subdivided in two thanks of a symmetry condition, this leads to a cells number of about 200 000. The minimum space mesh size is about 0.1 m.

Metas numerical setup

Concerning Metas, no steady state is performed before the train displacement. Therefore, the train-tunnel distance have to be long enough to obtain this steady state before the train entry. In addition, if the train-tunnel distance is too short, the sudden train displacement implies a contamination of the pressure signal. The initial train-tunnel distance is 160 m, against 15 m for Tg-flo. Thus, the domain is twice as long as the Tg-flo domain. The overall physical domain is modeled and no symmetry is used.

For the first calculation, meshes contain about 500 000 elements, and the minimum space mesh size is 0.2 m. This mesh is globally more coarse than Tg-flo mesh.

The calculation Metas2 was performed to obtain a better result. To realize that, Metas1 mesh is refined at two important locations: train nose and tunnel entry. With this refinement, the cells number reaches 1 360 000 and the minimum space mesh size is as for Tg-flo equal to 0.1 m.


Emmanuel Creusé
email: creuse@math.univ-lille1.fr

Serge Nicaise
email: serge.nicaise@univ-valenciennes.fr

A posteriori error estimator based on gradient recovery by averaging for discontinuous Galerkin methods

Keywords: A posteriori estimator, Discontinuous Galerkin finite elements AMS (MOS) subject classification 65N30; 65N15, 65N50

We consider some (anisotropic and piecewise constant) diffusion problems in domains of R 2 , approximated by a discontinuous Galerkin method with polynomials of any fixed degree. We propose an a posteriori error estimator based on gradient recovery by averaging. It is shown that this estimator gives rise to an upper bound where the constant is one up to some additional terms that guarantee reliability. The lower bound is also established. Moreover these additional terms are negligible when the recovered gradient is superconvergent. The reliability and efficiency of the proposed estimator is confirmed by some numerical tests.

Introduction

Among other methods, the finite element method is one of the more popular that is commonly used in the numerical realization of different problems appearing in engineering applications, like the Laplace equation, the Lamé system, the Stokes system, the Maxwell system, etc.... (see [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Ciarlet | The finite element method for elliptic problems[END_REF][START_REF] Monk | A posteriori error indicators for Maxwell's equations[END_REF]). More recently discontinuous Galerkin finite element methods become very attractive since they present some advantages. For example, they allow to do some "p refinement", by locally increasing the polynomial degree of the approximation if needed. They can moreover use non-conform meshes allowing hanging-nodes, making the mesh generation easier for concrete industrial cases. In our days a vast literature exists on the subject, we refer to [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF][START_REF] Cockburn | The development of discontinuous Galerkin methods[END_REF] and the references cited there. Adaptive techniques based on a posteriori error estimators have become indispensable tools for such methods. For continuous Galerkin finite element methods, there now exists a vast amount of literature on a posteriori error estimation for problems in mechanics or electromagnetism and obtaining locally defined a posteriori error estimates. We refer to the monographs [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] Babuška | The finite element methods and its reliability[END_REF][START_REF] Monk | Finite element methods for Maxwell's equations[END_REF][START_REF] Verfurth | A review of a posteriori error estimation and adaptive mesh-refinement techniques[END_REF] for a good overview on this topic. On the other hand a similar theory for discontinuous methods is less developed, let us quote [START_REF] Becker | Energy norm a posteriori error estimation for discontinuous Galerkin methods[END_REF][START_REF] Creusé | Anisotropic a posteriori error estimation for the mixed discontinuous Galerkin approximation of the Stokes problem[END_REF][START_REF] Houston | Energy norm a posteriori error estimation for mixed discontinuous Galerkin approximations of the Maxwell operator[END_REF][START_REF] Houston | A posteriori error estimation for discontinuous Galerkin discretization of the H(curl)-ellipic partial differential operator[END_REF][START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order problems[END_REF][START_REF] Rivière | A posteriori error estimates for a discontinuous Galerkin method applied to elliptic problems[END_REF][START_REF] Sun | L 2 (H 1 ) norm a posteriori error estimation for discontinuous Galerkin approximations of reactive transport problems[END_REF].

Usually upper and lower bounds are proved in order to guarantee the reliability and the efficiency of the proposed estimator. Most of the existing approaches involve constants depending on the shape regularity of the elements and/or of the jumps in the coefficients; but these dependences are often not given. Only a small number of approaches gives rise to estimates with explicit constants, let us quote [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF][START_REF] Ladevèze | Error estimate procedure in the finite element method and applications[END_REF][START_REF] Luce | A local a posteriori error estimator based on equilibrated fluxes[END_REF][START_REF] Neittaanmaäki | Reliable methods for computer simulation: error control and a posteriori error estimates[END_REF][START_REF] Nicaise | An a posteriori error estimator for the Lamé equation based on equilibrated fluxes[END_REF][START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF] for continuous methods. For discontinuous methods, we may cite the recent papers [START_REF] Ainsworth | A posteriori error estimation for discontinuous Galerkin finite element approximation[END_REF][START_REF] Lazarov | Functional a posteriori error estimates for discontinuous galerkin approximations of elliptic problems[END_REF][START_REF] Cochez-Dhondt | Equilibrated error estimators for discontinuous Galerkin methods[END_REF][START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF][START_REF] Ern | Guaranteed and robust discontinuous galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF].

Our goal is therefore to consider second order elliptic operators with discontinuous diffusion coefficients in two-dimensional domains with mixed boundary conditions and a discontinuous Galerkin method with polynomials of any degree. Inspired from the paper [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF], which treats the case of continuous diffusion coefficients approximated by a continuous Galerkin method, we further derive an a posteriori estimator with an explicit constant in the upper bound (more precisely 1) up to some additional terms that are usually superconvergent and some oscillating terms. The approach, called gradient recovery by averaging [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF] is based on the construction of a Zienkiewicz/Zhu estimator, namely the difference in an appropriate norm of a∇ h u h -Gu h , where ∇ h u h is the broken gradient of u h and Gu h is a H(div )-conforming approximation of this variable. Here special attention has to be paid due to the assumption that a may be discontinuous. Moreover the non conforming part of the error is managed using a Helmholtz decomposition of the error and a standard Oswald interpolation operator [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order problems[END_REF][START_REF] Ainsworth | A posteriori error estimation for discontinuous Galerkin finite element approximation[END_REF]. Furthermore using standard inverse inequalities, we show that our estimator is locally efficient. Two interests of this approach are first the simplicity of the construction of Gu h , and secondly its superconvergence property (validated by numerical tests).

The schedule of the paper is as follows: We recall in section 2 the diffusion problem, its numerical approximation and an appropriate Helmholtz decomposition of the error. Section 3 is devoted to the introduction of the estimator based on gradient averaging and the proofs of the upper and lower bounds. The upper bound directly follows from the construction of the estimator and some results from [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF], while the lower bound requires the use of some inverse inequalities and a special construction of Gu h . Finally in section 4 some numerical tests are presented that confirm the reliability and efficiency of our estimator and the superconvergence of Gu h to a∇u.

Let us finish this introduction with some notation used in the remainder of the paper: On D, the L 2 (D)-norm will be denoted by • D . In the case D = Ω, we will drop the index Ω. The usual norm and semi-norm of H s (D) (s ≥ 0) are denoted by • s,D and | • | s,D , respectively. Finally, the notation a b and a ∼ b means the existence of positive constants C 1 and C 2 , which are independent of the mesh size and of the quantities a and b under consideration such that a C 2 b and C 1 b a C 2 b, respectively. In other words, the constants may depend on the aspect ratio of the mesh, the diffusion matrix, as well as the polynomial degree l (see below).

The boundary value problem and its discretization

Let Ω be a bounded open domain of R 2 with a Lipschitz boundary Γ that we suppose to be polygonal. We further assume that Ω is simply connected and that Γ is connected. We consider the following elliptic second order boundary value problem with non homogeneous mixed boundary conditions:

   -div (a ∇u) = f in Ω, u = g D on Γ D , a∇u • n = g N on Γ N , (1)
where Γ = ΓD ∪ ΓN and Γ D ∩ Γ N = ∅. For convenience we suppose that meas Γ D > 0.

In the sequel, we suppose that a is a symmetric positive definite matrix which is piecewise constant, namely we assume that there exists a partition P of Ω into a finite set of Lipschitz polygonal domains Ω 1 , • • • , Ω J such that, on each Ω j , a = a j where a j is a symmetric positive definite matrix. The variational formulation of (1) involves the bilinear form

B (u, v) = Ω a∇u • ∇v
and the Hilbert space

H 1 D (Ω) = {u ∈ H 1 (Ω) : u = 0 on Γ D }. Given f ∈ L 2 (Ω), g D ∈ H 1 2 (Γ D ) and g N ∈ L 2 (Γ N ), the weak formulation consists in finding u ∈ w + H 1 D (Ω) such that B (u, v) = Ω f v + Γ N g N v, ∀v ∈ H 1 D (Ω), (2)
where w ∈ H 1 (Ω) is a lifting for g D , i.e., w = g D on Γ D . Invoking the positiveness of a, the bilinear form B is coercive on H 1 D (Ω) with respect to the norm

Ω |a 1/2 ∇u| 2 1/2
and this coerciveness guarantees that problem (2) has a unique solution by the Lax-Milgram lemma.

Discontinuous Galerkin approximated problem

Following [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order problems[END_REF][START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF], we consider the following discontinuous Galerkin approximation of our continuous problem: We consider a triangulation T of Ω, that is a "partition" of Ω made of triangles T (closed subsets of Ω) whose edges are denoted by e. We assume that this triangulation is regular, i.e., for any element T , the ratio h T ρ T is bounded by a constant σ > 0 independent of T and of the mesh size h = max T ∈T h T , where h T is the diameter of T and ρ T the diameter of its largest inscribed ball. We further assume that T is conforming with the partition P of Ω, i.e., the matrix a being constant on each T ∈ T , we then denote by a T the value of a restricted to an element T . With each edge e of the triangulation, we denote by h e its length, we associate a unit normal vector n e (whose orientation can be arbitrary chosen) and the so-called patch ω e = ∪ e∈T T , the union of triangles having e as edge. For a triangle T , n T stands for the outer unit normal vector of T . E (resp. N ) represents the set of edges (resp. vertices) of the triangulation. In the sequel, we need to distinguish between edges included into Ω, Γ D or Γ N , in other words, we set

E int = {e ∈ E : e ⊂ Ω}, E D = {e ∈ E : e ⊂ Γ D }, E N = {e ∈ E : e ⊂ Γ N }.
For shortness, we also write 2) is approximated by the (discontinuous) finite element space:

E ID = E int ∪ E D . Problem (
X h = v h ∈ L 2 (Ω)|v h|T ∈ P ℓ (T ), T ∈ T , (3)
where ℓ is a fixed positive integer.

For our further analysis we need to define some jumps and means of elements of X h through any e ∈ E of the triangulation. For e ∈ E such that e ⊂ Ω, we denote by T + and T -the two elements of T containing e. Let q ∈ X h , we denote by q ± , the traces of q taken from T ± , respectively. Then we define the mean of q on e by {{q}} = q + + q - 2 .

For v ∈ [X h ] d , we denote similarly

{{v}} = v + + v - 2 .
The jump of q on e is now defined as follows:

[[q]] = q + n T + + q -n T -.
Remark that the jump [[q]] of q is vector-valued. For a boundary edge e, i. e., e ⊂ ∂Ω, there exists a unique element T + ∈ T such that e ⊂ ∂T + . Therefore the mean and jump of q are defined by {{q}} = q + and [[q]] = q + n T + .

For q ∈ X h , we define its broken gradient ∇ h q in Ω by :

(∇ h q) |T = ∇q |T , ∀T ∈ T .
The space X h is equipped with the norm

q DG,h := T ∈T a 1/2 ∇q 2 T 1/2 + e∈E ID h -1 e [[q]] 2 e 1/2
.

Later on we also need the continuous counterpart of X h , namely we introduce

S h = v h ∈ C(Ω)|v h|T ∈ P ℓ (T ), T ∈ T ,
as well as

S h,1 = v h ∈ C(Ω)|v h|T ∈ P 1 (T ), T ∈ T .
We further need

X h,1 = v h ∈ L 2 (Ω)|v h|T ∈ P 1 (T ), T ∈ T .
With these notations, we define the bilinear form B h (., .) as follows:

B h (u h , v h ) := T ∈T T a∇u h • ∇v h - e∈E ID e ({{a∇ h v h }} • [[u h ]] + {{a∇ h u h }} • [[v h ]]) + e∈E ID γ e h -1 e e [[u h ]] • [[v h ]], ∀u h , v h ∈ X h ,
where the positive parameters γ e are chosen large enough to ensure coerciveness of the bilinear form B h on X h (see, e.g., Lemma 2.1 of [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order problems[END_REF]), namely according to Theorem 3 of [START_REF] Warburton | On the constants in hp-finite element trace inverse inequalities[END_REF] (and the arguments from section 3 of [START_REF] Shahbazi | An explicit expression for the penalty parameter of the interior penalty method[END_REF]), the following choices yield the coerciveness of B h :

γ e > c(ℓ + 1)(ℓ + 2) max T ⊂ωe C T , (4)
where C T is the largest eigenvalue of the diffusion matrix a T and c is the positive constant such that h 2 e ≤ c|T |, ∀T ⊂ ω e , that exists due to the regularity assumption on the mesh.

The discontinuous Galerkin approximation of problem ( 2) reads now: Find

u h ∈ X h , such that B h (u h , v h ) = F (v h ), (5)
where

F (v h ) = Ω f v h + e∈E D e g D (γh -1 e v h -a∇v h • n T ) + Γ N g N v h , ∀v h ∈ X h .
As our approximated scheme is a non conforming one (i.e. the solution does not belong to H 1 D (Ω)), as usual we need to use an appropriate Helmholtz decomposition of the error (see Lemma 3.2 of [START_REF] Dari | A posteriori error estimators for nonconforming finite element methods[END_REF] or Theorem 1 of [START_REF] Ainsworth | A posteriori error estimation for discontinuous Galerkin finite element approximation[END_REF] in 2D).

Lemma 2.1 (Helmholtz decomposition of the error) We have the following error decomposition

a∇ h (u -u h ) = a∇ϕ + curl χ, (6)
with χ ∈ H 1 (Ω) such that curl χ • n = 0 on Γ N , (7)
and ϕ ∈ H 1 D (Ω). Moreover the next identity holds:

a 1/2 ∇ h (u -u h ) 2 = a 1/2 ∇ h ϕ 2 + a -1/2 curl χ 2 . ( 8
)
For the detailed proof of this Lemma we refer to Lemma 2.1 of [START_REF] Cochez-Dhondt | Equilibrated error estimators for discontinuous Galerkin methods[END_REF].

The a posteriori error analysis based on gradient recovery by averaging

Error estimators can be constructed in many different ways as, for example, using residual type error estimators which measure locally the jump of the discrete flux [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order problems[END_REF]. A different method, based on equilibrated fluxes, consists in solving local Neumann boundary value problems [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF] or in using Raviart-Thomas interpolant [START_REF] Ainsworth | A posteriori error estimation for discontinuous Galerkin finite element approximation[END_REF][START_REF] Cochez-Dhondt | Equilibrated error estimators for discontinuous Galerkin methods[END_REF][START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF][START_REF] Ern | Guaranteed and robust discontinuous galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF]. Here, as an alternative we introduce a gradient recovery by averaging and define an error estimator based on a H(div )-conforming approximation of this variable. In comparison with [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF], we here allow the case of discontinuous diffusion coefficient and use a discontinuous Galerkin method.

Inspired from [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF] the conforming part of the estimator η CF involves the difference between the broken gradient a∇ h u h and its smoothed version Gu h , where Gu h is for the moment any element in X 2 h,1 satisfying

Gu h ∈ H(div , Ω) = {v ∈ L 2 (Ω) 2 : div v ∈ L 2 (Ω)}, ( 9
) (Gu h )| |Ω j ∈ H 1 (Ω j ), ∀j = 1, • • • , J. (10)
Hence conforming part of the estimator η CF is defined by

η 2 CF = T ∈T η 2 CF,T , (11)
where the indicator η CF,T is defined by

η CF,T = a -1/2 (a∇u h -Gu h ) T .
For the nonconforming part of the error, we associate with u h , its Oswald interpolation operator, namely the unique element w h ∈ S h defined in the following natural way (see Theorem 2.2 of [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order problems[END_REF]): to each node n of the mesh corresponding to the Lagrangian-type degrees of freedom of S h , the value of w h is the average of the values of u h at this node

n if it belongs to Ω ∪ Γ N (i.e., w h (n) = n∈T |T |u h|T (n) n∈T |T |
) and the value of g D at this node if it belongs to ΓD (here we assume that g D ∈ C( ΓD )). Then the non conforming indicator η N C,T is simply

η N C,T = a 1/2 ∇(w h -u h ) T .
The non conforming part of the estimator is then

η 2 N C = T ∈T η 2 N C,T . (12)
Similarly we introduce the estimator corresponding to jumps of u h :

η 2 J = e∈E ID η 2 J,e , with η 2 J,e = 1 he [[u h ]] 2 e if e ∈ E int , 1 he u h -g D 2 e
if e ∈ E D .

As in [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF], we introduce some additional superconvergent security parts. In order to define them properly we recall that for a node x ∈ N , we denote by λ x the standard hat function (defined as the unique element in S h,1 such that λ x (y) = δ x,y for all y ∈ N ), let ω x be the patch associated with x, which is simply the support of λ x and let h x be the diameter of ω x (which is equivalent to the diameter h K of any triangle K included into ω x ). We now denote by r the element residual

r = f + div (Gu h )
and for all x ∈ N , we set

rx = ( ωx λ x ) -1 ωx rλ x if x ∈ N \ N D , rx = 0 if x ∈ N D .
We further use a multilevel decomposition of S h,1 , namely we suppose that we start from a coarse grid T 0 and that the successive triangulations are obtained by using the bisection method, see [START_REF] Schmidt | Design of adaptive finite element software[END_REF][START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF]. This means that we obtain a finite sequence of nested triangulations

T ℓ , ℓ = 0, • • • , L such that T L = T . Denoting by S ℓ the space S ℓ = v ∈ C(Ω)|v |T ∈ P 1 (T ), T ∈ T ℓ ,
then we have S ℓ ⊂ S ℓ+1 and S h,1 = ∪ L ℓ=0 S ℓ = S L . Furthermore if we denote by N ℓ the nodes of the triangulation T ℓ , we have

N ℓ ⊂ N ℓ+1 .
As usual for all z ∈ N ℓ we denote by λ ℓz the hat function associated with z, namely the unique element in S ℓ such that

λ ℓz (z ′ ) = δ zz ′ ∀z ′ ∈ N ℓ .
For all ℓ ≥ 1 we finally set

Ñℓ = (N ℓ \ N ℓ-1 ) ∪ {z ∈ N ℓ-1 : λ ℓz = λ ℓ-1z },
and Ñ0 = N 0 . It should be noticed (see for instance [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF]) that to each z ∈ Ñℓ , the corresponding hat function λ ℓz does not belong to S ℓ-1 .

Now we define ρ and γ by

ρ2 = x∈N ρ 2 x , γ2 = L ℓ=0 z∈ Ñℓ \Γ D γ 2 ℓz ,
where

ρ 2 x = h 2 x ωx |r -rx | 2 λ x + h x ωx∩Γ N |Gu h • n -g N | 2 λ x , γ ℓz = | R, λ ℓz |,
R being the residual defined by

R, ϕ = Ω Gu h • ∇ϕ - Ω f ϕ - Γ N g N ϕ, ∀ϕ ∈ H 1 (Ω).

Upper bound

Theorem 3.1 Assume that there exists

w h ∈ X h ∩ H 1 (Ω) such that g D = w h|Γ D . Let u ∈ w + H 1 D (Ω)
be a solution of problem (2) and let u h be its discontinuous Galerkin approximation, i.e. u h ∈ X h solution of [START_REF] Becker | Energy norm a posteriori error estimation for discontinuous Galerkin methods[END_REF]. Then there exists C > 0 such that

a 1/2 ∇ h (u -u h ) ≤ (η 2 CF + η 2 N C ) 1/2 + C(ρ + γ), (13)
and consequently

u -u h DG,h ≤ (η 2 CF + η 2 N C ) 1/2 + η J + C(ρ + γ). ( 14
)
Remark 3.2 Let us note that under a superconvergence property of ||a -1/2 (Gu h -a∇u)||, ρ and γ will be proved to be negligible quantities (see Theorem 3.5 below), so that the error is in this case asymptotically bounded above by the estimator without any multiplicative constant. This superconvergence property is observed in most of practical cases, as for example in our numerical tests (see section 4). Moreover, theoretical results for different continuous finite element methods on structured and unstructured meshes have been established (see for example [START_REF] Hoffmann | Asymptotically exact a posteriori estimators for the pointwise gradient error on each element in irregular meshes, Part I : A smooth problem and globally quasi-uniform meshes[END_REF][START_REF] Schatz | Asymptotically exact a posteriori estimators for the pointwise gradient error on each element in irregular meshes, Part II : The piecewise linear case[END_REF][START_REF] Zhang | A posteriori error estimates on irregular grids based on gradient recovery[END_REF]), but, to our knowledge, not yet for discontinuous methods on unstructured multi-dimensional meshes.

Proof: From the Helmholtz decomposition of the error we have

a 1/2 ∇ h (u -u h ) 2 = a 1/2 ∇ϕ 2 + a -1/2 curl χ 2 . ( 15
)
We are then reduced to estimate each term of this right-hand side.

For the non conforming part, we proceed as in [START_REF] Ainsworth | A posteriori error estimation for discontinuous Galerkin finite element approximation[END_REF], namely by Green's formula we have

a -1/2 curl χ 2 = Ω ∇ h (u -u h ) • curl χ = - Ω ∇ h u h • curl χ + Γ D g D curl χ • n = Ω ∇ h (w h -u h ) • curl χ, since Ω ∇w h • curl χ = Γ D g D curl χ • n.
By Cauchy-Schwarz's inequality we directly obtain

a -1/2 curl χ 2 ≤ η N C a -1/2 curl χ . ( 16
)
For the conforming part, we write

a 1/2 ∇ϕ 2 = Ω a∇ h (u -u h ) • ∇ϕ = Ω (a∇u -Gu h ) • ∇ϕ + Ω (Gu h -a∇ h u h ) • ∇ϕ.
By Cauchy-Schwarz's inequality we obtain

a 1/2 ∇ϕ 2 ≤ a -1/2 (a∇ h u h -Gu h ) a 1/2 ∇ϕ + Ω (a∇u -Gu h ) • ∇ϕ . ( 17
)
Using problem (2), the second term of this right-hand side is bounded by the residual, indeed

- Ω (a∇u -Gu h ) • ∇ϕ = Ω Gu h • ∇ϕ - Ω f ϕ - Γ N g N ϕ = R, ϕ .
Using the arguments from Theorem 4.1 of [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF], we have

| R, ϕ | ≤ C(ρ + γ) a 1/2 ∇ϕ . ( 18
)
Coming back to the identity [START_REF] Ern | Guaranteed and robust discontinuous galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF], and using the estimates ( 16), ( 17) and ( 18) we conclude by discrete Cauchy-Schwarz's inequality and again using (15):

a 1/2 ∇ h (u -u h ) 2 ≤ η N C a -1/2 curl χ + (η CF + C(ρ + γ)) a 1/2 ∇ϕ ≤ (η 2 N C + η 2 CF ) 1/2 ( a -1/2 curl χ 2 + a 1/2 ∇ϕ 2 ) 1/2 + C(ρ + γ) a 1/2 ∇ϕ ≤ [(η 2 N C + η 2 CF ) 1/2 + C(ρ + γ)] a 1/2 ∇ h (u -u h ) .

Lower bound

Our lower bound is based on the equivalence of the local L 2 -norm of any element in X h with a local L 2 -norm in the interfaces. First of all for any vertex x of one Ω j and belonging to more than one sub-domain, we introduce the following local notation: let Ω i , i = 1, • • • , n, n ≥ 2, the sub-domains that have x as vertex. We further denote by n i the unit normal vector along the interface I i between Ω i and Ω i+1 (modulo n if x is inside the domain Ω) and oriented from Ω i and Ω i+1 . Now we are able to prove the following lemma:

Lemma 3.3 Assume that x is a vertex of one Ω j and belonging to more than one subdomain, and use the notations introduced above. Then there exists a positive constant C that depends only on the geometrical situation of the Ω i 's near x such that for all

v (i) ∈ R 2 , i = 1, • • • , n, there exist vectors g(v) (i) ∈ R 2 , i = 1, • • • , n satisfying (g(v) (i+1) -g(v) (i) ) • n i = 0, ∀i = 1, • • • , n, (19)
and such that the following estimate holds

n i=1 |v (i) -g(v) (i) | ≤ C n i=1 |[[v • n]] i |, (20)
where here | • | means the Euclidean norm and

[[v • n]] i means the normal jump of v along the interface I i : [[v • n]] i = (v (i+1) -v (i) ) • n i , ∀i = 1, • • • , n.
Proof: First introduce the following subspace of R 2n :

W = {v = (v (i) ) n i=1 : v (i) ∈ R 2 and satisfying [[v • n]] i = 0, ∀i = 1, • • • , n}.
We take g(v) = Π W v, the orthogonal projection of v = (v (i) ) n i=1 into W . By construction g(v) trivially satisfies [START_REF] Hoffmann | Asymptotically exact a posteriori estimators for the pointwise gradient error on each element in irregular meshes, Part I : A smooth problem and globally quasi-uniform meshes[END_REF]. On the other hand the estimate ( 20) is equivalent to

|v -Π W v| ≤ C n i=1 |[[v • n]] i |,
which is easily proved by a contradiction argument and the fact that we are in a finite dimensional space.

Using the above lemma, we are now able to prove the asymptotic nondeterioration of the smoothed gradient if the following choice for Gu h is made: We distinguish the following different possibilities for x ∈ N . 1) First for all vertex x of the mesh (i.e. vertex of at least one triangle) such that x is inside one Ω j , we set

(Gu h ) |Ω j (x) = 1 |ω x | x∈T |T |a T ∇u h|T (x). (21)
2) Second if x belongs to the boundary of Ω and to the boundary of only one Ω j (hence it does not belong to the boundary of another Ω k ), we define (Gu h ) |Ω j (x) as before.

3) If x belongs to an interface between two different sub-domain Ω j and Ω k but is not a vertex of these sub-domains, then we denote by n j,k the unit normal vector pointing from Ω j to Ω k and set t j,k the unit orthogonal vector of n j,k so that (n j,k , t j,k ) is a direct basis of R 2 ; in that case we set

(Gu h ) |Ω j (x) • n j,k = (Gu h ) |Ω k (x) • n j,k = 1 |ω x | x∈T |T |a T ∇u h|T (x) • n j,k , (22)
(Gu h ) |Ω j (x) • t j,k = 1 |ω x ∩ Ω j | T ⊂Ω j :x∈T |T |a T ∇u h|T (x) • t j,k , (23)
(Gu h ) |Ω k (x) • t j,k = 1 |ω x ∩ Ω k | T ⊂Ω k :x∈T |T |a T ∇u h|T (x) • t j,k . (24)
4) Finally if x is a vertex of at least two sub-domains Ω j , for the sake of simplicity we suppose that each triangle T having x as vertex is included into one Ω j , and we take

(Gu h ) |Ω j (x) = g(v) (j) ∀j ∈ J x , (25)
where

J x = {j ∈ {1, • • • , J} : x ∈ Ωj }, g(v) (j)
were defined in the previous Lemma 3.3 with here v given by v = (a j ∇u h|T (x)) j∈Jx . With these choices, we take

(Gu h ) |Ω j = x∈N ∩ Ωj (Gu h ) |Ω j (x)λ x , ∀j = 1, • • • , J, (26)
where (Gu h ) |Ω j (x) was defined before.

The main point is that by construction Gu h satisfies the requirements ( 9) and ( 10) but moreover we have the next asymptotic nondeterioration result: 3)), then for each element T ∈ T the following estimate holds

Theorem 3.4 If ℓ ≤ 2 (see (
a -1/2 T (Gu h -a T ∇u) T u -u h DG,ω T + osc(f, ω T ), ( 27
)
where ω T denotes the patch consisting of all the triangles of T having a nonempty intersection with T and

v 2 DG,ω T = a 1/2 ∇ h v 2 ω T + γ e∈E ID :e⊂ω T h -1 e [[v]] 2 e , and
osc(f, ω T ) 2 = h 2 T T ′ ⊂ω T f -Π T ′ f 2 T ′ ,
where Π T ′ f is the L 2 (T ′ )-orthogonal projection of f onto P 1 (T ′ ).

Proof: By the triangle inequality we may write

a -1/2 T (Gu h -a T ∇u) T ≤ a -1/2 T (Gu h -a T ∇u h ) T + a -1/2 T (a T ∇u h -a T ∇u) T ≤ a -1/2 T (Gu h -a T ∇u h ) T + u -u h DG,T .
Therefore it remains to estimate the first term of this right-hand side. For that purpose, since T ⊂ Ω j for a unique j ∈ {1, • • • , J}, we may write

(Gu h -a T ∇u h ) |T = x∈T {(Gu h ) |Ω j (x) -a j ∇u h|T (x)}λ x .
As 0 ≤ λ x ≤ 1, and since the triangulation is regular, we get

a -1/2 T (Gu h -a T ∇u h ) T x∈T |(Gu h ) |Ω j (x) -a j ∇u h|T (x)|h T . ( 28
)
We are then reduced to estimate the factor |(Gu h ) |Ω j (x)a j ∇u h|T (x)| for all nodes x of T . For that purpose, we distinguish four different cases: 1) If x ∈ Ω j , then we use an argument similar to the one from Proposition 4.2 of [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF] adapted to the DG situation. By the definition of Gu h , we have

(Gu h ) |Ω j (x) = 1 |ω x | T ′ ⊂ωx |T ′ |a j ∇u h|T ′ (x),
because in this case all T ′ ⊂ ω x are included into Ω j . As a consequence, we obtain

(Gu h ) |Ω j (x) -a j ∇u h|T (x) = 1 |ω x | T ′ ⊂ωx |T ′ |a j (∇u h|T ′ (x) -∇u h|T (x)),
and therefore

|(Gu h ) |Ω j (x) -a j ∇u h|T (x)| T ′ ⊂ωx |a j (∇u h|T ′ (x) -∇u h|T (x))|.
For each T ′ ⊂ ω x , there exists a path of triangles of ω x , written T i , i = 0, • • • , n such that

T 0 = T, T n = T ′ , T i = T j , ∀i = j, T i ∩ T i+1 is an common edge ∀i = 1, • • • , n -1.
Hence by the triangle inequality we can estimate

|a j (∇u h|T ′ (x) -∇u h|T (x))| ≤ n-1 i=0 |a j (∇u h|T i+1 (x) -∇u h|T i (x))|.
Now for each term, since a j is symmetric and positive definite, we have

|a j (∇u h|T i+1 (x)-∇u h|T i (x))| |{a j (∇u h|T i+1 (x)-∇u h|T i (x))}•n i |+|(∇u h|T i+1 (x)-∇u h|T i (x))•t i |,
where n i is one fixed unit normal vector along the edge T i ∩ T i+1 and t i is one fixed unit tangent vector along this edge.

All together we have shown that

|(Gu h ) |Ω j (x) -a j ∇u h|T (x)|h T h T e∈E int :e⊂ωx {|[[a j ∇u h (x) • n]] e | + |[[∇ -hu h (x) • t]] e |}.
Using a norm equivalence and an inverse inequality we obtain

|(Gu h ) |Ω j (x) -a j ∇u h|T (x)|h T e∈E int :e⊂ωx {h 1/2 e [[a j ∇ h u h • n]] e e + h -1/2 e [[u h ]] e }. (29)
2) If the node x belongs to the boundary of Ω and to the boundary of a unique Ω j , since (Gu h ) |Ω j (x) is defined as in the first case, the above arguments lead to [START_REF] Schmidt | Design of adaptive finite element software[END_REF].

3) If x is a vertex of different sub-domains Ω j , then by Lemma 3.3, we have

|(Gu h ) |Ω j (x) -a j ∇u h|T (x)|h T h T e∈E int :e⊂ωx |[[a∇ h u h (x) • n]] e |, (30)
and therefore as before we conclude that (29) holds. 4) Finally if x belongs to an interface between two subdomains and is not a vertex of them, then it is not difficult to show that (30) holds (due to the regularity of the mesh), and consequently ( 29) is still valid. Summarizing the different cases, by ( 28) and ( 29), we have

a -1/2 T (Gu h -a T ∇u h ) T e∈E int :e⊂ωx {h 1/2 e [[a∇ h u h • n]] e e + h -1/2 e [[u h ]] e }. (31)
The first term of this right hand side is a part of the standard residual error estimator and it is by now standard that (using appropriate bubble functions and Green's formula)

h 1/2 e [[a∇u h • n]] e e a∇ h (u -u h ) ωe + osc(f, ω e ), ∀e ∈ E int .
The second term is part of the DG-norm. Therefore the above estimate in (31) leads to [START_REF] Nicaise | An a posteriori error estimator for the Lamé equation based on equilibrated fluxes[END_REF]. Now using the same arguments than in Proposition 4.1 of [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF], we have Theorem 3.5 For all T ∈ T , x ∈ N and ℓ ≥ 0, z ∈ N ℓ , we have

η CF,T ≤ a 1/2 T ∇(u h -u) T + a -1/2 (Gu h -a∇u) T , ρx a -1/2 (Gu h -a∇u) ωx + osc(f, ω x ) + osc(g N , ω x ), γℓz a -1/2 (Gu h -a∇u) ω ℓz ,
where

osc(g N , ω x ) 2 = e⊂ωx∩Γ N h e g N -Π e g N 2
T
, Π e g N being the L 2 (e)-orthogonal projection of g N onto P 2 (e).

For the non conforming part of the estimator, we make use of Theorem 2.2 of [START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order problems[END_REF] to directly obtain the Theorem 3.6 Let the assumptions of Theorem 3.1 be satisfied. For each element T ∈ T the following estimate holds

η N C,T a 1/2 T u -u h DG,ω T . (32)
A direct consequence of these three Theorems is the next local lower bound:

Theorem 3.7 Let the assumptions of Theorems 3.1 and 3.4 be satisfied. For each element T ∈ T the following estimate holds

η CF,T + η N C;T + η J,T + x∈T (ρ x + γx ) u -u h DG,ω T + osc(f, ω T ) + osc(g N , ω T ),
where γx = γLx recalling that L is such that N L = N .

Remark 3.8 Note that the lower bound on the non conforming estimator (see [START_REF] Schatz | Asymptotically exact a posteriori estimators for the pointwise gradient error on each element in irregular meshes, Part II : The piecewise linear case[END_REF]) involves a constant that depends on the aspect ratio of the mesh and of the penalization parameter γ, and is specific to the discontinuous Galerkin method. Consequently it prevents the estimator to be asymptotically exact, as in the continuous Galerkin method [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF].

Nevertheless the numerical tests show quite satisfactory effectivity indices (see below).

As in Proposition 4.3 of [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF], one has γ a -1/2 (Gu h -a∇u) ,

and therefore a global lower bound can be obtained: Theorem 3.9 Let the assumptions of Theorems 3.1 and 3.4 be satisfied. Then the following global lower bound holds

η CF + η N C + η J + ρ + γ u -u h DG,h + osc(f, Ω) + osc(g N , Ω).
4 Numerical results

The polynomial solution

In order to illustrate our theoretical predictions, this first numerical test consists in validating our computations on a simple case, using an uniform refinement process. Let Ω be the square (-1, 1) 2 , Γ D = ∂Ω, a = Id and f defined such that the exact solution u is given by : u(x, y) = (x + 1)(x -1)(y + 1)(y -1).

Let us recall that u h is the finite element solution, e(u h ) = ||uu h || DG,h the error,

η(u h ) = η 2 CF + η 2 N C
1/2 + η J the estimator and Gu h the approximated value of a∇u given by [START_REF] Neittaanmaäki | Reliable methods for computer simulation: error control and a posteriori error estimates[END_REF]. We also define CV error (resp. CV recov ) as the convergence rate of the error e(u h ) (resp. of the quantity ||Gu h -∇u||) with respect to DoF -1/2 from one line of the table to the following one.

Computations are performed using a global mesh refinement process from an initial cartesian grid, using γ e ≡ γ = 20 to ensure (4). First, it can be seen from Table 1 that the convergence rate of the numerical method is equal to one, as theoretically expected. Then, the superconvergence property of the term ||Gu h -∇u|| is actually observed with a convergence rate of 1.50. Finally, the reliability of the estimator is ensured since the ratio in the last column (the so-called effectivity index), converges fastly towards the constant 1.70.

k

The interior and boundary layer solution

The following numerical test consists in solving the interior and boundary layer example given in [START_REF] Fierro | A posteriori error estimators, gradient recovery by averaging, and superconvergence[END_REF]. Let Ω the square (-1, 1) 2 , Γ D = ∂Ω, a = Id and f defined such that :

u(x, y) = arctan(60(x 2 + y 2 -1))
is the exact solution (see Figure 1). Let us note that the boundary layer crosses the boundary ∂Ω and that the loading term oscillates across it, what constitutes the difficulty of the computation. This time, an adaptive mesh refinement strategy is used based on the estimator η T = η CF,T + η N C,T + η J,T and the marking procedure

η T > 0.75 max T ′ η T ′
and a standard refinement procedure with a limitation on the minimal angle. Once again, we choose γ e ≡ γ = 20. Several mesh levels are displayed on Figure 2, to show the capability of the algorithm to track the high gradients regions. Furthermore, quantitative results are displayed on Table 2. Once again, the superconvergence property of ||Gu h -∇u|| is asymptotically observed, as well as the reliability of the estimator, provided the mesh is fine enough around the boundary layer.

The discontinuous case

This section is devoted to the case of the discontinuous coefficient a. Namely, the domain Ω = (-1, 1) 2 with Γ D = Γ is decomposed into 4 sub-domains Ω i , i = 1, ..., 4, with Ω 1 = (0, 1) × (0, 1), Ω 2 = (-1, 0) × (0, 1), Ω 3 = (-1, 0) × (-1, 0) and Ω 4 = (0, 1) × (-1, 0). In that case we take discontinuous coefficient a, namely we take a = a i on Ω i , with a 1 = a 3 = Id and a 2 = a 4 = C Id with C to be specified. For this second test, and using usual polar coordinates centered at (0, 0), the exact solution is equal to the singular function u(x, y) = r α φ(θ), where α ∈ (0, 1) and φ are chosen such that u is harmonic on each sub-domain Ω i , i = 1, ..4, and satisfies the jump conditions :

[u] = 0 and [a∇u.n] = 0

h,1 = v h ∈ C(Ω) 2 |v h|T ∈ P 2 1 (T ), T ∈ T and characterized by its value at each node of the mesh given by :

(G * u h )(x) = 1 |ω x | x∈T |T |a T ∇u h|T (x). ( 33
)
For the case C = 5, we choose γ e ≡ γ = 20 and figure 3 shows some of the meshes obtained during the local refinement process, which is the same as the one used in section 4.2. Moreover, table 3 displays the corresponding quantitative results. The smoothed gradient Gu h is superconvergent, while the effectivity index converges towards 1.75, which is quite satisfactory and comparable with results from [START_REF] Cochez-Dhondt | Equilibrated error estimators for discontinuous Galerkin methods[END_REF][START_REF] Ern | Guaranteed and robust discontinuous galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF] as well as those of the previous tests in sections 4.1 and 4.2.

For the case C = 100, we choose γ e ≡ γ = 500 and figure 4 shows some of the meshes obtained during the local refinement process. These ones are more refined around the singularity than the case C = 5. Table 4 displays the corresponding quantitative results. Once again, the smoothed gradient Gu h is superconvergent, provided the mesh resolution is fine enough. The effectivity index requires more iterations to converge, but is already equal to 2.82 for the finer mesh resolution.

Let us finally note that for both cases, the standard ZZ smoothed gradient G * u h is no more superconvergent towards a∇u, as in the case a = Id. This is not surprising since the statement of Theorem 3.4 is not valid for G * u h if a = Id. k k

Figure 1 :

 1 Figure 1: Interior and boundary layer solution.

Figure 2 :

 2 Figure 2: Mesh levels 1, 10 and 20, interior and boundary layer solution.

Figure 3 :

 3 Figure 3: Mesh levels 1, 3 and 10, singular solution for C = 5.

Figure 4 :

 4 Figure 4: Mesh levels 1, 3 and 10, singular solution for C = 100.

Table 1 :

 1 The polynomial solution (uniform refinement).

	DoF	e(u h )	CV error ||Gu h -∇u|| CV recov	η(u h ) e(u h )
	1 384	3.80E-01		2.78E-01		1.60
	2 1536	1.90E-01	1.00	1.04E-01	1.42	1.67
	3 6144	9.51E-02	1.00	3.73E-02	1.47	1.70
	4 24576 4.74E-02	1.00	1.32E-02	1.49	1.70
	5 98304 2.37E-02	1.00	4.69E-03	1.50	1.70
	6 393216 1.19E-02	1.00	1.66E-02	1.50	1.70

Table 2 :

 2 The interior and boundary layer solution (adaptive refinement).

	k DoF	e(u h )	CV error ||Gu h -∇u|| CV recov	η(u h ) e(u h )
	1 384	1.02E+02			4.67E+01		1.10
	7 1323	4.59E+01	1.29	3.36E+01	0.53	0.94
	15 6069	1.71E+01	1.29	1.67E+01	0.91	1.35
	22 22458 7.84E+00	1.19	6.55E+00	1.43	1.67
	28 94260 3.56E+00	1.10	1.84E+00	1.77	1.62
	34 366180 1.80E+00	1.01	7.14E-01	1.40	1.60
	on the interfaces. Non-homogeneous Dirichlet boundary conditions on Γ are fixed accord-
	ingly. It is easy to see (see for instance [11]) that α is the root of the transcendental
	equation		tan	α π 4	= 1/C.

Since α < 1, this solution has a singular behavior around the point (0, 0). For this test, we also compute the standard ZZ smoothed gradient G * u h belonging to S

Table 3 :

 3 DoF e(u h )CVerror ||Gu h -∇u|| CV recov ||G * u h -∇u|| e(u h) Discontinuous coefficient a: C = 5, γ = 20, singular solution (local refinement).

	η(u h )
	e(u h )

Table 4 :

 4 DoF e(u h )CV error ||Gu h -∇u|| CV recov ||G * u h -∇u|| e(u h ) Discontinuous coefficient a: C = 100, γ = 500, singular solution (local refinement).

	η(u h )
	e(u h )