Sarah Delcourte

Nathalie Glinsky

ANALYSIS OF A DISCONTINUOUS GALERKIN METHOD FOR ELASTODYNAMIC EQUATIONS. APPLICATION TO 3D WAVE PROPAGATION

Keywords: AMS subject classifications: 35L50, 35F10, 35F15, 35L05, 35Q99 discontinuous Galerkin method, centered flux, leap-frog scheme, elastodynamic equation

In this paper, we introduce a second-order leap-frog time scheme combined with a high-order discontinuous Galerkin method for the solution of the 3D elastodynamic equations. We prove that this explicit scheme is stable under a CFL type condition obtained from a discrete energy which is preserved in domains with free surface or decreasing in domains with absorbing boundary conditions. Moreover, we study the convergence of the method for both the semi-discrete and the fully discrete scheme, and we illustrate the convergence results by the propagation of an eigenmode. Finally, we examine a more realistic 3D test case simulating the propagation of the wave produced by an explosive source in a half-space which constitutes a validation of the source introduction and the absorbing boundary conditions.

Introduction

In the last few decades, the computational seismology has become an essential tool to simulate realistic wavefields of local, regional or even global wave propagation problems. The physics governing these phenomena is now well understood and many different accurate numerical methods have been developed and can deal with threedimensional realistic applications thanks to a continuous increase of the computational resources and the use of parallel computational facilities.

Among all the numerical methods proposed for simulations in time domain, the most popular is undoubtedly the finite difference (FD) method and its many improvements from the initial FD schemes proposed by Alterman and Karal [START_REF] Alterman | Propagation of elastic waves in layered media by finite-difference methods[END_REF] or Kelly et al. [START_REF] Kelly | Synthetic seismograms : a finite-difference approach[END_REF] such as, for instance, the introduction of the velocity-stress system and the staggeredgrids (Madariaga [START_REF] Madariaga | Dynamics of an expanding circular fault[END_REF], Virieux [START_REF] Virieux | P-SV wave propagation in heterogeneous media, velocity-stress finite difference method[END_REF]), fourth-order schemes in space (Bayliss et al. [START_REF] Bayliss | A fourth-order accurate finite-difference scheme for the computation of elastic waves[END_REF], Levander [START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF]) and rotated staggered-grids (Saenger et al. [START_REF] Saenger | Modeling the propagation of elastic waves using a modified finite-difference grid[END_REF]) allowing strong fluctuations of the elastic parameters. Their major drawback is the restriction to cartesian grids not suited for geometrical internal or surfacic nonlinearities (topography) but, the use of non uniform grids (Pitarka [START_REF] Pitarka | 3D elastic finite-difference modeling of seismic motion using staggered-grids with nonuniform spacing[END_REF]) or even discontinuous grids (Aoi and Fujiwara [START_REF] Aoi | 3D Finite-difference method using discontinuous grids[END_REF]) enabled improving the accuracy at the free surface and, more recently, an hybrid method coupling a finite-difference technique for the most part of the domain and a finite-element method in subdomains containing the nonlinearities (topography, faults) has been proposed (Moczo et al. [START_REF] Moczo | Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures[END_REF], Galis et al. [START_REF] Galis | A 3-D hybrid finite-difference-finite-element viscoelastic modelling of seismic wave motion[END_REF]). Some other methods have been further developed such as finite element (FE) methods which allow meshes adapted to complex geometries [START_REF] Lysmer | A finite element method for seismology[END_REF], [START_REF] Marfurt | Accuracy of finite-difference and finite-element modelling of the scalar and elastic wave equations[END_REF], [START_REF] Bao | Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers[END_REF]. However, they are very costly because one needs to invert a global mass matrix at each time step. This difficulty was overcome by the use of Gauss-Lobato Legendre quadrature formulae at the root of the spectral element methods (SEM) which have been widely applied to quadrangular and hexaedral meshes (see, for instance [START_REF] Komatitsch | The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures[END_REF][START_REF] Chaljub | Vilotte Solving elastodynamics in a fluid-solid heterogeneous sphere: a parallel spectral element approximation on non-conforming grids[END_REF] amongst many contributions). Note the recent development of two non conforming discontinuous Galerkin spectral element methods and their convergence and stability analysis for the wave propagation in 2D quadrangular meshes [START_REF] Antonietti | Non-conforming high order approximations of the elastodynamics equation[END_REF] and 3D affine hexahedral meshes [START_REF] Bui-Thanh | Analysis of an hp-nonconforming discontinuous Galerkin spectral element method for wave propagation[END_REF]. In very complex media, the use of simplicial meshes permits a better approximation of the geometry (bann, topography, faults). For this reason, we study a high-order Discontinuous Galerkin (DG) method applied to tetrahedral meshes.

The Discontinuous Galerkin (DG) method has been initially introduced by Reed and Hill [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF] for the solution of neutron transport problems. Neglected during twenty years, it became very popular to solve hyperbolic problems especially computational electromagnetics [START_REF] Fezoui | Convergence and stability of a Discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF]. In spite of its success in many domains of applications, this method has been rarely applied to seismic wave propagation problems. Käeser et al. (see [START_REF] Käser | Quantitative accuracy analysis of the discontinuous Galerkin method for seismic wave propagation[END_REF][START_REF] Dumbser | An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes II: the three-dimensional isotropic case[END_REF] and references therein) proposed a DG finite element scheme that uses the ADER approach based on an upwind scheme in order to solve the elastic wave equations with the same high accuracy in both space and time. However, despite its accuracy, this scheme is diffusive; it is the reason why we propose another combination which preserves a discrete elastodynamic energy.

In this paper, we study the P-SV seismic wave propagation considering an isotropic, linearly elastic medium by solving the velocity-stress formulation of the elastodynamic equations. For the discretization of this system, we focus on a DG method which is a finite element method allowing discontinuities at the interfaces introduced via numerical fluxes as for finite volumes. Our method is based on centered fluxes and a leap-frog time-discretization which lead to a non-dissipative combination [START_REF] Delcourte | A high-order Discontinuous Galerkin method for the seismic wave propagation[END_REF]. Moreover, the method is suitable for complex unstructured simplicial meshes. The extension to higher order in space is realized by Lagrange polynomial functions (of degree 0 to 3 for our solver), defined locally on tetrahedra and do not necessitate the inversion of a global mass matrix since an explicit time scheme is used.

This article is organized as follows. In section 2, we state the velocity-stress formulation in a symmetrical pseudoconservative form. Then, in section 3, we detail the discretization of the equations system by the Discontinuous Galerkin method based on centered fluxes in space, combined with a leap-frog scheme in time. The approximation of the boundary conditions is also presented. After, we study, in section 4, the preservation of a discrete elastodynamic energy and the stability of the scheme, taking into account the free surface or the absorbing boundary conditions. Section 5 is devoted to the convergence analysis of the semi-discrete and fully discretized schemes. Finally, we illustrate this study in section 6 by some numerical results in three dimensions of space: the propagation of an eigenmode and the case of an explosive source in a half-space.

Velocity-stress formulation in pseudo-conservative form

The P-SV wave propagation in an isotropic, linearly elastic medium is modelised by the elastodynamic equations, which initially write in displacement-stress formulation; let be U = (U α ) α∈{x,y,z} the displacement vector and σ = (σ α,β ) α,β∈{x,y,z} , the stress tensor, then the system reads

   ρ ∂ 2 t U = ∇ • σ , σ = λ (∇ • U ) I + µ (∇ U + (∇ U ) t ) , (1)
where I is the identity matrix, ρ is the density of the medium and λ and µ are the Lamé parameters related to shear and compressional velocities (v s and v p ) in the medium by v s = µ ρ and v p = λ+2µ ρ . We introduce the velocity vector V = ∂ t U in the equation ( 1) and obtain the first-order velocity-stress formulation [START_REF] Virieux | P-SV wave propagation in heterogeneous media, velocity-stress finite difference method[END_REF]:

   ρ ∂ t V = ∇ • σ , ∂ t σ = λ (∇ • V ) I + µ (∇ V + (∇ V ) t ) . (2)
Let W = V , σ t be the vector composed of the velocity components V = (V x , V y , V z ) t and the stress components σ = (σ xx , σ yy , σ zz , σ xy , σ xz , σ yz ) t , since the stress tensor is symmetrical, then, the system (2) can be rewritten as

∂ t W + α∈{x,y,z} A α (ρ, λ, µ) ∂ α W = 0 . (3)
We choose here not to detail the matrices A α (ρ, λ, µ).

In order to express the system (3) in a pseudo conservative form, we introduce the following change of variables on the stress components,

σ = T σ = σ xx + σ yy + σ zz 3 , 2σ xx -σ yy -σ zz 3 , -σ xx + 2σ yy -σ zz 3 , σ xy , σ xz , σ yz . (4)
So, the system reads

       ρ ∂ t V = α∈{x,y,z} M α ∂ α σ , Λ 0 (λ, µ) ∂ t σ = α∈{x,y,z} N α ∂ α V , (5)
where

Λ 0 (λ, µ) = diag 3 3λ + 2µ , 3 2µ , 3 2µ , 1 µ , 1 µ , 1
µ is a diagonal matrix containing the characteristics of the medium and

M x =   1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0   , M y =   0 0 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 1   , M z =   0 0 0 0 1 0 0 0 0 0 0 1 1 -1 -1 0 0 0   , N x =         1 0 0 2 0 0 -1 0 0 0 1 0 0 0 1 0 0 0         , N y =         0 1 0 0 -1 0 0 2 0 1 0 0 0 0 0 0 0 1         , N z =         0 0 1 0 0 -1 0 0 -1 0 0 0 1 0 0 0 1 0         .
We can notice that now the matrices M α and N α (α = x, y, z) are constant and do not depend anymore on the material properties. So, the system ( 5) is a pseudo-conservative formulation of (3). At last, we multiply the second equation of ( 5) by the following positive definite symmetrical (PDS) matrix

S =         1 0 0 0 0 0 0 2 3 1 3 0 0 0 0 1 3 2 3 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1        
in order to obtain a symmetrical system. Therefore, we finally get the symmetrical pseudo-conservative formulation:

       ρ∂ t V = α∈{x,y,z} M α ∂ α σ , Λ(λ, µ)∂ t σ = α∈{x,y,z} M t α ∂ α V , (6)
where Λ(λ, µ) = SΛ 0 (λ, µ) is a positive definite symmetrical matrix whose spectrum is

Sp{Λ} = 1 2µ , 3 2µ , 1 µ , 1 µ , 1 µ , 3 3λ + 2µ ⊂ R * +
and we can notice that SN α = M t α (α = x, y, z). This formulation will be very useful to establish the energy preservation. Then, we add: (i) a physical boundary condition on the free surface

σ n = 0, (7)
which will be rewritten by respect to σ in section 3.3.

(ii) absorbing boundary conditions

P σ = -A V and P t V = -B σ, (8)
(where the matrices P, A and B will be specified in section 3.3) to approximate an infinite domain.

At last, we need to impose initial data

V (0) = V 0 and σ(0) = σ0 , (9)
to ensure the existence and the uniqueness of the solution ( V , σ) of the system ( 6)-( 7)-( 8)-( 9).

3 A discontinuous Galerkin method combined with a leap-frog scheme

Integration on a simplex

We consider a bounded polyhedral domain Ω of R 3 , discretized in N T tetrahedra T i , which form a partition of the domain. We assume that the characteristics of the medium are constant over each element T i and denoted by (ρ i , λ i , µ i ). We multiply, in the sense of the scalar product, the first (resp. the second) equation of the problem (6) by a vector field ϕ i ∈ R 3 (resp. ψ i ∈ R 6 ) and we integrate them on each element T i . Then, we apply the Green formula to the right-hand sides:

ρ i Ti ϕ t i ∂ t V dv = α∈{x,y,z} - Ti (∂ α ϕ i ) t M α σdv + ∂Ti ϕ t i M α σn αi ds , ( 10a
) Ti ψ t i Λ i ∂ t σdv = α∈{x,y,z} - Ti (∂ α ψ i ) t M t α V dv + ∂Ti ψ t i M t α V n αi ds , (10b)
where we set Λ i = Λ(λ i , µ i ) and n i = (n xi , n yi , n zi ) t represents the outwards unit normal vector to T i .

Evaluation of volume integrals

The approximations of V and σ are denoted by the fields V h and σh which are defined locally on each element T i and may be discontinuous through the interfaces, so that we set

∀i, V h | T i = V i and σh | T i = σi . ( 11
)
Consider the Lagrange nodal interpolants ϕ ij , ψ ij ∈ P m (T i ), set of polynomials over T i with a degree m, then we can write V i and σi as linear combinations of time-dependent fields:

V i (x, y, z, t) = dof j=1 V ij (t) ϕ ij (x, y, z) , (12)
σi (x, y, z, t)

= dof j=1 σij (t) ψ ij (x, y, z) , (13)
where dof denotes the number of degrees of freedom on the element T i . To approximate the volume integrals on T i , we just have to replace the fields V and σ in (10) by V h and σh .

Approximation on faces

To calculate the integrals on ∂T i of (10), we split this boundary in internal and boundary faces. We denote by V(i) the set of indices of the neighboring elements of T i and we note by S ik each internal face common to both elements T i and T k , i.e. S ik = T i ∩ T k . Moreover, some elements T i have one or more faces common to the boundary of the domain. The set of the indices k of such faces S ik = T i ∩ ∂Ω is denoted by E(i) for absorbing boundaries and K(i)

for free surface boundaries. Remark that, for most elements, E(i) and K(i) are empty sets.

ψ t i M t α V n αi ds = k∈V(i) S ik ψ t i P t ik V ds + k∈E(i) S ik ψ t i P t ik V ds + k∈K(i) S ik ψ t i P t ik V ds , (14b)
where n ik = (n x ik , n y ik , n z ik ) t represents the unit normal vector of S ik , oriented from T i towards T k .

For the interior faces S ik of T i (k ∈ V(i)), we choose to apply centered fluxes by introducing the mean-value on this face:

V |S ik ≃ 1 2 V i + V k and σ|S ik ≃ 1 2 σi + σk . ( 15
)
For the boundary integrals, two types of boundary conditions have been considered : a free surface condition at the physical interface between air and the medium, and an absorbing condition on the artificial boundaries of an infinite domain.

Free surfaces: On these faces, we compute the fluxes by introducing weakly the physical condition [START_REF] Benjemaa | Etude de stabilité d'un schéma volumes finis pour les équations de l'élasto-dynamique en maillages non structurés[END_REF] in the third term of (14a) and the physical condition is rewritten in the new variable σ via the change of variables [START_REF] Bao | Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers[END_REF].

No physical condition is applied to the velocity in the third term of (14b). So, for a boundary face S ik of T i , with k ∈ K(i), this condition reduces to:

V |S ik ≃ V i and P ik σ|S ik = 0 . (16)
Absorbing surfaces: To simulate infinite domains, we introduce artificial boundaries, on which we impose absorbing conditions. Therefore, for any k ∈ E(i) and any real unit vector n ik , we define the matrix

A n ik (ρ i , λ i , µ i ) = α∈{x,y,z} A α (ρ i , λ i , µ i ) n α ik = -   0 R 3×3 1 ρ i P ik Λ -1 i P t ik 0 R 6×6   , (17)
which is diagonalizable in R, i.e. all its eigenvalues e k (k=1,...,9) are real:

e 1 = -v pi ,
v pi = v p (ρ i , λ i , µ i ) = λ i + 2µ i ρ i and v si = v s (ρ i , λ i , µ i ) = µ i ρ i (18)
and we note by P n (ρ, λ, µ) the matrix whose column k is the right eigenvector associated to the eigenvalue e k :

0 B B B B B B B B B B B B B B B B B B B B B B @ 3vpnx -vsnxny vsn0,1,-1 0 0 0 -vsn0,1,-1 vsnxny -3vpnx 3vpny vsn1,0,-1 -vsnxny 0 0 0 vsnxny -vsn1,0,-1 -3vpny 3vpnz vsnynz vsnxnz 0 0 0 -vsnxnz -vsnynz -3vpnz 3λ + 2µ 0 0 -nxn0,1,1 -nyn1,0,1 -nzn1,1,0 0 0 3λ + 2µ 2µn2,-1,-1 -2µn 2 x ny 2µnxn0,1,-1 nxn0,1,1 nyn1,0,-2 nzn1,-2,0 2µnxn0,1,-1 -2µn 2 x ny 2µn2,-1,-1 -2µn1,-2,1 2µnyn1,0,-1 -2µnxn 2 y nxn0,1,-2 nyn1,0,1 nzn-2,1,0 -2µnxn 2 y 2µnyn1,0,-1 -2µn1,-2,1 6µnxny µnxn1,-1,-1 µny n-1,1,-1 0 0 3nxnynz µnyn-1,1,-1 µnxn1,-1,-1 6µnxny 6µnxnz 0 µnzn1,1,-1 0 3nxnynz 0 µnz n1,1,-1 0 6µnxnz 6µnynz µnz n1,1,-1 0 3nxnynz 0 0 0 µnz n1,1,-1 6µnynz 1 C C C C C C C C C C C C C C C C C C C C C C A
where, to simplify, we set n a,b,c = an 2 x + bn 2 y + cn 2 z . Then, the boundary absorbing conditions consist in an upwind technique where we only take into account the outgoing waves neglecting the ingoing waves part. For that, we approximate A n ik by the matrix A + n ik where we set A + n ik = P n ik D + P -1 n ik (ρ i , λ i , µ i ) with D + (ρ i , λ i , µ i ) the diagonal matrix composed by the positive eigenvalues of A n ik . Therefore, according to [START_REF] Glinsky | A high-order discontinuous Galerkin scheme for elastic wave propagation[END_REF], we have to approximate

P ik σ P t ik V = - ρ i I 3 0 R 3×6 0 R 6×3 Λ i A n ik (ρ i , λ i , µ i ) V σ , (19)
where I 3 denotes the identity matrix, by

P ik σ P t ik V ≃ - ρ i I 3 0 R 3×6 0 R 6×3 Λ i A + n ik (ρ i , λ i , µ i ) V i σi . (20)
After calculus, we obtain:

ρ i I 3 0 R 3×6 0 R 6×3 Λ i A + n ik (ρ i , λ i , µ i ) = 1 2 A ik -P ik -P t ik B ik (21)
where

A ik = ρ i   (v pi -v si )n 2 x ik + v si (v pi -v si )n x ik n y ik (v pi -v si )n x ik n z ik (v pi -v si )n x ik n y ik (v pi -v si )n 2 y ik + v si (v pi -v si )n y ik n z ik (v pi -v si )n x ik n z ik (v pi -v si )n y ik n z ik (v pi -v si )n 2 z ik + v si   (22)
and

B ik = 1 ρ i v si         1 n 2 x ik -n 2 z ik n 2 y ik -n 2 z ik 2n x ik n y ik 2n x ik n z ik 2n y ik n z ik n 2 x ik -n 2 z ik n 2 x ik + n 2 z ik n 2 z ik n x ik n y ik 0 -n y ik n z ik n 2 y ik -n 2 z ik n 2 z ik n 2 y ik + n 2 z ik n x ik n y ik -n x ik n z ik 0 2n x ik n y ik n x ik n y ik n x ik n y ik n 2 x ik + n 2 y ik n y ik n z ik n x ik n z ik 2n x ik n z ik 0 -n x ik n z ik n y ik n z ik n 2 x ik + n 2 z ik n x ik n y ik 2n y ik n z ik -n y ik n z ik 0 n x ik n z ik n x ik n y ik n 2 y ik + n 2 z ik         + v si -v pi ρ i v pi v si         1 n 2 x ik -n 2 z ik n 2 y ik -n 2 z ik 2n x ik n y ik 2n x ik n z ik 2n y ik n z ik         ⊗         1 n 2 x ik -n 2 z ik n 2 y ik -n 2 z ik 2n x ik n y ik 2n x ik n z ik 2n y ik n z ik         . (23)
Finally, for a boundary face S ik of T i with k ∈ E(i), the fluxes are approximated in [START_REF] Galis | A 3-D hybrid finite-difference-finite-element viscoelastic modelling of seismic wave motion[END_REF] by:

S ik ϕ t i P ik σ ds ≃ S ik ϕ t i 2 -A ik V i + P ik σi ds , (24)
S ik

ψ t i P t ik V ds ≃ S ik ψ t i 2 P t ik V i -B ik σi ds . ( 25
)
Remark 1 The absorbing boundary conditions [START_REF] Benjemaa | Etude et simulation numérique de la rupture dynamique des séismes par des méthodes d'éléments finis discontinus[END_REF] follow from ( 20)- [START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF]. It is a first-order approximation, efficient for waves with a normal incidence to the artificial boundaries but whose accuracy could be improved especially to limit reflections in presence of grazing waves.

Remark 2 A ik is a symmetrical positive definite matrix whose spectrum is

Sp(A ik ) = {ρ i v pi , ρ i v si , ρ i v si } ,
whereas B ik is a symmetrical semi-definite positive matrix. Indeed, we have

det(B ik -xI 6 ) = x 3 (ax 3 + bx 2 + cx + d) where a = 1, b = - 1 ρ i v pi v si v pi (2 -2n 2 y ik + n 2 z ik -n 4 z ik + 2n 4 y ik + 2n 2 y ik n 2 z ik ) + v si (2 + n 4 z ik -2n 4 y ik + 2n 2 y ik -2n 2 y ik n 2 z ik ) , c = 1 ρ 2 i v pi v 2 si v pi (1 + 2n 4 y ik -2n 2 y ik -n 4 z ik + n 2 z ik -n 2 y ik n 2 z ik + 3n 2 y ik n 4 z ik + 3n 4 y ik n 2 z ik ) +v si (4 -n 2 y ik + n 4 y ik + 2n 2 y ik + 4n 2 y ik n 2 z ik -3n 2 y ik n 4 z ik -3n 4 y ik n 2 z ik + 3n 4 y ik n 2 z ik ) , d = - 1 ρ 3 i v pi v 2 si 2 + 2n 2 z ik -3n 2 y ik + 3n 4 y ik -n 4 z ik -n 2 y ik n 2 z ik + 4n 4 y ik n 2 z ik + 4n 2 y ik n 4 z ik .
We can deduce that 0 is a third order eigenvalue of B ik , but the three other eigenvalues are not easily specified. Then, we shall use the Sylvester's criteria. For that, we denote by B ik = (b jl ) 1≤j,l≤6 the coefficients of the matrix B ik and B q = (b jl ) 1≤j,l≤q the submatrices of B ik , and we observe that

det(B 1 ) = 1 ρ i v pi > 0, det(B 2 ) = 1 ρ 2 i v pi v si 4n 2 x ik n 2 z ik + n 2 y ik (n 2 x ik + n 2 z ik ) ≥ 0, det(B 3 ) = 9 ρ 3 i v pi v 2 si n 2 x ik n 2 y ik n 2 z ik ≥ 0, det(B 4 ) = det(B 5 ) = det(B 6 ) = 0,
so that the three other eigenvalues of B ik are necessarily positive.

Time-discretization

At last, we apply a leap-frog time-integration scheme. It is a two steps explicit scheme which results, when combined to the centered fluxes defined at [START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF], in a non-diffusive scheme (see section 4.2). We note by ∆t the time-step and the superscripts refer to time stations. Thus, on each element T i , V n i represents the velocity field at t n = n∆t and σn+ 1 2 i the stress tensor components at t n+ 1 2 = (n + 1 2 )∆t. Starting from initial values V 0 i at t = 0 and σ 1 2 i at t = ∆t 2 , the final scheme expresses on each element T i by:

ρ i Ti ϕ t i V n+1 i -V n i ∆t dv = - α∈{x,y,z} Ti (∂ α ϕ i ) t M α σn+ 1 2 i dv + k∈V(i) S ik ϕ t i P ik σn+ 1 2 i + σn+ 1 2 k 2 ds + k∈E(i) S ik ϕ t i 1 2 P ik σn+ 1 2 i - 1 2 A ik V n i ds , ( 26
) Ti ψ t i Λ i σn+ 3 2 i - σn+ 1 2 i ∆t dv = - α∈{x,y,z} Ti (∂ α ψ i ) t M t α V n+1 i dv + k∈V(i) S ik ψ t i P t ik V n+1 i + V n+1 k 2 ds + k∈E(i) S ik ψ t i 1 2 P t ik V n+1 i - 1 2 B ik σn+ 1 2 i ds + k∈K(i) S ik ψ t i P t ik V n+1 i ds . (27)
Obviously, linear combinations ( 12) and ( 13) are still valid with superscripts.

Energy preservation and stability of the scheme

The discrete scheme ( 26)-( 27) being explicit, it is conditionally stable. The aim of this section is to establish a sufficient condition on the time step ∆t for the L 2 -stability of the Discontinuous Galerkin scheme, taking into account both free surface ( 16) and absorbing boundary conditions ( 24)- [START_REF] Moczo | Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures[END_REF]. Following [START_REF] Fezoui | Convergence and stability of a Discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF][START_REF] Benjemaa | Etude de stabilité d'un schéma volumes finis pour les équations de l'élasto-dynamique en maillages non structurés[END_REF][START_REF] Benjemaa | Etude et simulation numérique de la rupture dynamique des séismes par des méthodes d'éléments finis discontinus[END_REF], we shall define a discrete energy and prove that it is a quadratic positive definite form playing the role of a Lyapunov function of all the unknowns V n i and σn-1/2 i .

A discrete elastodynamic energy

In the continuum, the total energy of the system is given by:

E(t) = 1 2 Ω ρ V (t) 2 dv E k , kinetic energy + 1 2 Ω σ(t) : ǫ(t) dv Em, mechanical energy
, where ǫ is the symmetrical deformation tensor whose components are ǫ i,j = 1 2 ∂Ui ∂xj + ∂Uj ∂xi . As for the stress, we introduce the deformation vector ǫ = (ǫ xx , ǫ yy , ǫ zz , 2 ǫ xy , 2 ǫ xz , 2 ǫ yz )

t . For an elastic medium, the generalized Hooke's law links deformations and stresses through the linear relationship ǫ = C σ, where the matrix C expresses

C = 1 E         1 -ν -ν 0 0 0 -ν 1 -ν 0 0 0 -ν -ν 1 0 0 0 0 0 0 2 (1 + ν) 0 0 0 0 0 0 2 (1 + ν) 0 0 0 0 0 0 2 (1 + ν)        
, E and ν being respectively the Young's modulus and the Poisson's ratio

E = µ (3λ + 2µ) λ + µ and ν = λ 2 (λ + µ) .
Thanks to this relationship, the mechanical energy becomes

E m = 1 2 Ω ( σ) t ǫ dv = 1 2 Ω ( σ) t C σ dv ,
and introducing the change of variables (4), we obtain

E m = 1 2 Ω T -1 σ t C T -1 σ dv = 1 2 Ω σ t T -1 t C T -1 σ dv .
We can easily check that T -1 t C T -1 = Λ so that, the total energy writes

E(t) = 1 2 Ω ρ V (t) 2 dv + 1 2 Ω σ(t) t Λ(λ, µ) σ(t) dv .
Definition 1 The discrete energy E n in three dimensions of space at time n∆t is defined by:

E n = 1 2 NT i=1 Ti ρ i V n i t V n i + σn+ 1 2 i t Λ i σn-1 2 i dv . ( 28
)

Energy preservation

Now, we study the evolution of the discrete elastodynamic energy through one time-step. More precisely, we aim to establish that the combination between the centered approximation of the fluxes and the leap-frog time-integration leads to a non-dissipative scheme for unbounded domains or domains with free surfaces.

For the simplex T i , we denote by

V [n+ 1 2 ] i = V n+1 i + V n i 2 and σ[n] i = σn+ 1 2 i + σn-1 2 i 2
the mean-values of the velocity and the stress vectors at times n + 1 2 ∆t and n∆t respectively. Lemma 1 The variation of the total discrete energy during one time-step ∆t is

E n+1 -E n = - ∆t 2 N T X i=1 X k∈E(i) Z S ik " " V [n+ 1 2 ] i « t A ik V n i + " σn+ 1 2 i « t B ik σ[n] i # ds . ( 29
)
Proof 1 We calculate the variation of the discrete energy during one time-step ∆t:

E n+1 -E n = NT i=1 Ti   ρi V n+1 i + V n i t 2 V n+1 i -V n i + σn+ 1 2 i t 2 Λ i σn+ 3 2 i - σn-1 2 i    dv . (30)
Thus, we substitute ϕ i and ψ i in ( 26)-( 27) by

V [n+ 1 2 ] i and σn+ 1 2 i
whose each component belongs to P m (T i ):

E n+1 -E n = -∆t NT i=1 α∈{x,y,z} Ti ∂ α V [n+ 1 2 ] i t M α σn+ 1 2 i + ∂ α σn+ 1 2 i t M t α V [n+ 1 2 ] i dv + ∆t 2 NT i=1 k∈V(i) S ik V [n+ 1 2 ] i t P ik σn+ 1 2 k + σn+ 1 2 i t P t ik V [n+ 1 2 ] k ds - ∆t 2 NT i=1 k∈E(i) S ik V [n+ 1 2 ] i t A ik V n i + σn+ 1 2 i t B ik σ[n] i ds + ∆t 2 NT i=1   k∈V(i) S ik V [n+ 1 2 ] i t P ik σn+ 1 2 i + σn+ 1 2 i t P t ik V [n+ 1 2 ] i ds + k∈E(i) S ik V [n+ 1 2 ] i t P ik σn+ 1 2 i + σn+ 1 2 i t P t ik V [n+ 1 2 ] i ds + 2 k∈K(i) S ik σn+ 1 2 i t P t ik V [n+ 1 2 ] i ds   . (31)
By the Green formula, we obtain:

Ti ∂ α V [n+ 1 2 ] i t M α σn+ 1 2 i + ∂ α σn+ 1 2 i t M t α V [n+ 1 2 ] i dv = ∂Ti V [n+ 1 2 ] i t M α σn+ 1 2 i n αi ds . (32)
Moreover, we clearly have

σn+ 1 2 i t P t ik V [n+ 1 2 ] i = V [n+ 1 2 ] i t P ik σn+ 1 2 i
. Consequently, we deduce that the first and the fourth terms of E n+1 -E n vanish. Further, thanks to the orientation of n α ik from T i to T k , we remark that:

NT i=1 k∈V(i) S ik σn+ 1 2 i t P t ik V [n+ 1 2 ] k ds = - NT i=1 k∈V(i) S ik V [n+ 1 2 ] i t P ik σn+ 1 2 k ds ,
which implies that the second term is also equal to zero. Therefore, it remains only the third term which concludes this proof.

Corollary 1 Using the scheme ( 26)-( 27), for an infinite domain or a domain including free surface boundaries ( 16) only, the discrete elastodynamic energy is preserved through one time-step:

∀n ∈ N, E n+1 = E n .
Proof 2 In an infinite domain or a domain with free surface boundary conditions, E(i) is an empty set. Therefore, following (29), the discrete elastodynamic energy is preserved.

A corrected discrete elastodynamic energy for absorbing boundaries

When absorbing boundary conditions are applied at some faces of the domain (i.e. when E(i) is not empty), the discrete variation of the energy through one time-step ∆t is not necessarily negative according to Lemma 1, because of the time-asymmetry on both V i and σi in [START_REF] Saenger | Modeling the propagation of elastic waves using a modified finite-difference grid[END_REF]. In order to overcome this difficulty, we introduce correction terms in the definition of the discrete energy [START_REF] Scheid | Convergence of a Discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF] and prove that this corrected discrete elastodynamic energy is not increased through each time step.

Definition 2 We introduce a corrected discrete elastodynamic energy E n at time n∆t by:

E n = E n - ∆t 8 NT i=1 k∈E(i) S ik V n i t A ik V n i - σn-1 2 i t B ik σn-1 2 i ds. (33)
These correction terms only concern the absorbing faces and have no particular physical meaning. However, they match the loss of energy (which appears as a decreasing discrete energy) through these faces since we only consider outgoing waves.

Lemma 2 The sequel (E n ) n∈N is not increasing. More precisely, we have for all n ∈ N:

E n+1 -E n = - ∆t 2 NT i=1 k∈E(i) S ik V [n+ 1 2 ] i t A ik V [n+ 1 2 ] i + σ[n] i t B ik σ[n] i ds. ( 34
)
Proof 3 We calculate the variation of the total corrected discrete elastodynamic energy:

E n+1 -E n = - ∆t 4 NT i=1 k∈E(i) S ik V n i + V n+1 i t A ik V n i + σn+ 1 2 i t B ik σn-1 2 i + σn+ 1 2 i ds - ∆t 8 NT i=1 k∈E(i) S ik V n+1 i t A ik V n+1 i -V n i t A ik V n i ds + ∆t 8 NT i=1 k∈E(i) S ik σn+ 1 2 i t B ik σn+ 1 2 i - σn-1 2 i t B ik σn-1 2 i ds .
After a reorganization of the terms by using the symmetry of the matrices A ik and B ik , the last expression simplifies into (34), which is negative because A ik and B ik are respectively a positive definite symmetrical and a semi-positive definite symmetrical matrices.

Stability of the scheme

We aim at proving that the discrete corrected elastodynamic energy E n is a positive definite quadratic form of the unknowns V n i and σn-1/2 i under some stability condition on the time step ∆t. For that, we need to introduce another formulation of E n , independently of the unknowns σn+1/2 i .

Lemma 3

The discrete corrected energy E n , defined at (33), can be rewritten as

E n = 1 2 NT i=1 Ti ρ i V n i t V n i + σn-1 2 i t Λ i σn-1 2 i dv + ∆t 4 NT i=1 α∈{x,y,z} Ti ∂ α V n i t M α σn-1 2 i -∂ α σn-1 2 i t M t α V n i dv + ∆t 4 NT i=1 k∈V(i) S ik σn-1 2 i t P t ik V n k ds + ∆t 4 NT i=1 k∈K(i) S ik σn-1 2 i t P t ik V n i ds - ∆t 8 NT i=1 k∈E(i) S ik V n i t A ik V n i + σn-1 2 i t B ik σn-1 2 i ds. ( 35
)
Proof 4 At first, we consider Eq. ( 27) at the previous time step (as the velocity is at time n∆t) and we replace ψ i by σn-

1 2 i : Ti σn-1 2 i t Λ i σn+ 1 2 i dv = Ti σn-1 2 i t Λ i σn-1 2 i dv -∆t α∈{x,y,z} Ti ∂ α σn-1 2 i t M t α V n i dv +∆t k∈V(i) S ik σn-1 2 i t P t ik V n i + V n k 2 ds +∆t k∈E(i) S ik σn-1 2 i t 1 2 P t ik V n i - 1 2 B ik σn-1 2 i ds +∆t k∈K(i) S ik σn-1 2 i t P t ik V n i ds.
Then, we substitute the previous expansion in the definition of the discrete corrected energy (33) and we get:

E n = 1 2 NT i=1 Ti   ρ i V n i t V n i + σn-1 2 i t Λ i σn-1 2 i -∆t α∈{x,y,z} ∂ α σn-1 2 i t M t α V n i   dv + ∆t 4 NT i=1 k∈E(i) S ik σn-1 2 i t P t ik V n i - 1 2 V n i t A ik V n i - 1 2 σn-1 2 i t B ik σn-1 2 i ds + ∆t 4 NT i=1 k∈V(i) S ik σn-1 2 i t P t ik V n i + V n k ds + ∆t 2 NT i=1 k∈K(i) S ik σn-1 2 i t P t ik V n i ds.
At last, we conclude this proof by applying the Green formula:

k∈V(i) S ik σn-1 2 i t P t ik V n i ds + k∈E(i) S ik σn-1 2 i t P t ik V n i ds + k∈K(i) S ik σn-1 2 i t P t ik V n i ds = α∈{x,y,z} Ti ∂ α σn-1 2 i t M t α V n i + ∂ α V n i t M α σn-1 2 i dv.
In what follows, we denote by S i = k∈V(i)

S ik + k∈K(i) S ik + k∈E(i)
S ik = 0 the area of a tetrahedron T i . Furthermore, we assume that inside each finite element T i , there exist two positive dimensionless constants a i > 0 and b ik > 0 (for k ∈ V(i), k ∈ K(i), k ∈ E(i)) such that, for all polynomial function φ defined on T i :

∂ α φ L 2 (Ti) ≤ a i S i |T i | φ L 2 (Ti) , (36)
φ 2 L 2 (S ik ) ≤ b ik S ik |T i | φ 2 L 2 (Ti) , (37)
where φ L 2 (Ti) and φ L 2 (S ik ) mean the L 2 -norm of φ over T i and over the face S ik respectively. Note that the constants a i and b ik do not depend on the size of the finite element T i because they are invariant by any homothetic transformation. In fact, they only depend on the geometry of T i and on the shapes of φ (see [START_REF] Fezoui | Convergence and stability of a Discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF]).

These assumptions will be applied to V n i and σn-

1 2 i
which are linear combinations of Lagrange polynomial basis functions.

Theorem 1 Considering the scheme ( 26)-( 27), under assumptions (36) and (37), the corrected total discrete elastodynamic energy E n is non-increasing through iterations and is a positive definite quadratic form of all the unknowns V n i and σn-1/2 i . Consequently, the scheme is L 2 -stable under the CFL type condition:

∀i ∈ [1, N T ], ∀k ∈ V(i) ∪ K(i) ∪ E(i) : ∆t [14a i + b ik max(7, r + (A ik ), r + (B ik ))] < 4 |T i | S i min (ρ i , r -(Λ i )), ( 38
)
where we denote by r -

(Λ i ) = min 1 2µ i , 3 3λ i + 2µ i > 0 the smallest eigenvalue of Λ i , r + (A ik ) = ρ i v pi > 0 the
greatest eigenvalue of A ik and r + (B ik ) > 0 the greatest eigenvalue of B ik .

Remark 3 B ik is a real symmetrical semi-definite positive matrix by the Sylvester criteria and the eigenvalue 0 has a multiplicity equal to three. Consequently, B ik has three strictly positive bounded eigenvalues depending only on the material properties (whose expansions are not simple) and then, r + (B ik ) > 0.

Proof 5 We can bound down the first term of the right-hand side of (35) in the following way:

Ti ρ i V n i t V n i + σn-1 2 i t Λ i σn-1 2 i dv ≥ ρ i V n i 2 L 2 (Ti) 3 + r -(Λ i ) σn-1 2 i 2 L 2 (Ti) 6 . (39)
Then, we set M = r(MM t ) where r(MM t ) is the spectral radius of the matrix MM t . Thus, M x = M y = 2 and M z = 3, so that P ik ≤ α∈{x,y,z} M α = 7. Applying (36), we clearly have:

α∈{x,y,z} Ti ∂ α V n i t M α σn-1 2 i -∂ α σn-1 2 i t M t α V n i dv ≤ 14
a i S i |T i | V n i L 2 (Ti) 3 σn-1 2 i L 2 (Ti) 6 ≤ 7 a i S i |T i | V n i 2 L 2 (Ti) 3 + σn-1 2 i 2 L 2 (Ti) 6 . ( 40
)
On the other hand, using (37) for the internal faces, we also have:

S ik σn-1 2 i t P t ik V n k ds ≤ P ik 2 V n k 2 L 2 (S ik ) 3 + σn-1 2 i 2 L 2 (S ik ) 6 ≤ 7 2 b ki S ki |T k | V n k 2 L 2 (T k ) 3 + b ik S ik |T i | σn-1 2 i 2 L 2 (Ti) 6 . (41)
In the very same way, for the boundary faces with free surface condition, we have:

S ik σn-1 2 i t P t ik V n i ds ≤ 7 2 b ik S ik |T i | V n i 2 L 2 (Ti) 3 + σn-1 2 i 2 L 2 (Ti) 6 . (42)
Then, the last term of the right-hand side of (35) can be bounded thanks to:

S ik V n i t A ik V n i + σn-1 2 i t B ik σn-1 2 i ds ≤ r + (A ik ) V n i 2 L 2 (S ik ) 3 + r + (B ik ) σn-1 2 i 2 L 2 (S ik ) 6 ≤ b ik S ik |T i | r + (A ik ) V n i 2 L 2 (Ti) 3 + r + (B ik ) σn-1 2 i 2 L 2 (Ti) 6 .
Therefore, as S i = k∈V(i)

S ik + k∈K(i) S ik + k∈E(i)
S ik = 0 and thanks to the following remark related with (41):

NT i=1 k∈V(i) b ki S ki |T k | V n k 2 L 2 (T k ) 3 = NT i=1 1 |T i | V n i 2 L 2 (Ti) 3 k∈V(i) b ik S ik , (43)
we finally obtain the following estimation:

E n ≥ 1 2 NT i=1 V n i 2 L 2 (Ti) 3   k∈V(i) S ik ρ i S i - 7 4 ∆t |T i | (2a i + b ik ) + k∈K(i) S ik ρ i S i - 7 4 ∆t |T i | (2a i + b ik ) + k∈E(i) S ik ρ i S i - 1 4 ∆t |T i | (14a i + b ik r + (A ik ))   + 1 2 NT i=1 σn-1 2 i 2 L 2 (Ti) 6   k∈V(i) S ik r -(Λ i ) S i - 7 4 ∆t |T i | (2a i + b ik ) + k∈K(i) S ik r -(Λ i ) S i - 7 4 ∆t |T i | (2a i + b ik ) + k∈E(i) S ik r -(Λ i ) S i - 1 4 ∆t |T i | (14a i + b ik r + (B ik )   .
Finally, for any face S ik of the mesh, the discrete corrected energy is a positive definite quadratic form of all the unknowns V n i and σn-

1 2 i
under the condition (38) on the time step.

Remark 4 In the condition (38), the ratio

|T i | S i
has the same dimension as the diameter h i of the finite element T i , so that (38) is a CFL type stability condition.

Remark 5

The condition (38) might be suboptimal. Indeed, when studying this energy in the particular case of the finite volume P 0 approximation, for a uniform mesh with reference tetrahedra (in this case, a i = 0 and b ik = 1, ∀i ∈ [1, N T ]), in an infinite domain, we obtain a more precise CFL condition:

v p ∆t h < √ 2 3 √ 3
,

where h = max i∈[1,NT ]
h i and the medium is supposed to be homogeneous.

On the other hand, we remark that max(7, r + (A ik ), r + (B ik )) ≥ r + (A ik ) = ρv p and min (ρ i , r -(Λ i )) ≤ ρ, so that in some cases, (38) is simplified by:

∀i ∈ [1, N T ], v p ∆t < 4 |T i | S i , (44)
where

|T i | = h 3 √ 2 24
and

S i = h 2 (3 + √ 3) 4 .

Remark 6

The stability result (38) remains valid when the degree of the polynomial approximation is defined locally and may vary from an element T i to its neighbour T k .

Convergence analysis

The objective of this section is to prove the convergence of the totally discretized scheme ( 26)-( 27) following the method in [START_REF] Fezoui | Convergence and stability of a Discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes[END_REF][START_REF] Scheid | Convergence of a Discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF]. We consider a family of unstructured tetrahedral meshes T h , where h = max i∈ [START_REF] Antonietti | Non-conforming high order approximations of the elastodynamics equation[END_REF]NT ] h i , which forms a partition of the domain Ω, i.e. Ω = ∪ Ti∈T h T i . We assume that the unstructured meshes T h are uniformly shape regular in the sense that there exists a constant ξ > 0 such that

∀h > 0, ∀T i ∈ T h : h i d i ≤ ξ , (45)
where d i is the diameter of the biggest ball included in T i . Moreover, we assume that there exists η > 0 (independent of h) such that

∀h > 0, ∀T i ∈ T h , ∀k ∈ V(i) : h k h i ≤ η. (46)
In what follows, we denote by

X 1 h = { ϕ h ∈ L 2 (Ω) 3 : ∀i, ϕ h | T i ∈ P k (T i ) 3 } , (47)
X 2 h = { ψ h ∈ L 2 (Ω) 6 : ∀i, ψ h | T i ∈ P k (T i ) 6 } , (48)
the sets of discontinuous piecewise polynomial fields of degree at most k on each T i . We next introduce the broken Sobolev spaces

P H s (Ω) = v : ∀i, v |T i ∈ H s (T i ) (49) equipped with the norm v P H s (Ω) = NT i=1 v |T i 2 H s (Ti) 1/2
, where • H s (Ti) denotes the standard H s -norm on T i .

Definition and properties of the semi-discretized scheme

First, we are interested in the study of consistency and stability of the spatially semi-discretized problem below.

In what follows, we set ϕ i (resp.

ψ i ) the restriction of ϕ ∈ H s (T h ) 3 (resp. ψ ∈ H s (T h ) 6 ) to the element T i with s > 1 2 . The semi-discrete solution ( V h , σh ) defined in C 1 ([0, T ], X 1 h × X 2 h
) is the solution of the weak formulation: for any test field (

ϕ h , ψ h ) ∈ X 1 h × X 2 h , for 0 ≤ t ≤ T and all i, ρ i Ti ϕ t i ∂ t V i dv = - α∈{x,y,z} Ti (∂ α ϕ i ) t M α σi dv + k∈V(i) S ik ϕ t i P ik σi + σk 2 ds + k∈E(i) S ik ϕ t i 1 2 P ik σi - 1 2 A ik V i ds , (50)
Ti ψ t i Λ i ∂ t σi dv = - α∈{x,y,z} Ti (∂ α ψ i ) t M t α V i dv + k∈V(i) S ik ψ t i P t ik V i + V k 2 ds + k∈E(i) S ik ψ t i 1 2 P t ik V i - 1 2 B ik σi ds + k∈K(i) S ik ψ t i P t ik V i ds , (51)
with the initial values V h (0) = P 1 h ( V 0 ) and σh (0

) = P 2 h ( σ0 ), (52)
where P 1 h : L 2 (Ω) 3 → X 1 h and P 2 h : L 2 (Ω) 6 → X 2 h denote respectively the orthogonal projection onto X 1 h and X 2 h . Definition 3 For given vector fields W = ( V , σ) and W ′ = ( V ′ , σ′ ), we introduce the following bilinear forms which are well defined on

X 1 h × X 2 h : m( W , W ′ ) = NT i=1 Ti ρ i ( V ′ i ) t V i + ( σ′ i ) t Λ i σi dv (53) a( W , W ′ ) = NT i=1 Ti α∈{x,y,z} ∂ α V ′ i t M α σi + ∂ α σ′ i t M t α V i dv (54) b( W , W ′ ) = - NT i=1 k∈V(i) S ik ( V ′ i ) t P ik σi + σk 2 + ( σ′ i ) t P t ik V i + V k 2 ds - 1 2 NT i=1 k∈E(i) S ik ( V ′ i ) t P ik σi -A ik V i + ( σ′ i ) t P t ik V i -B ik σi ds - NT i=1 k∈K(i) S ik ( σ′ i ) t P t ik V i ds. (55)
Thanks to (50) and (51), we remark the following result:

Proposition 1 The semi-discrete solution W h = ( V h , σh ) of ( 50)-( 51)-(52) satisfies

m(∂ t W h , W ′ h ) + a( W h , W ′ h ) + b( W h , W ′ h ) = 0, ∀ W ′ h ∈ X 1 h × X 2 h . (56)
Proposition 2 The exact solution W = ( V , σ) of ( 6)-( 7)-( 8)-( 9) satisfies the following property

m(∂ t W , W ′ h ) + a( W , W ′ h ) + b( W , W ′ h ) = 0, ∀ W ′ h ∈ X 1 h × X 2 h . ( 57
)
Proof 6 According to (10a)-(10b), we have to verify that the quantity b( W , W ′ h ) is equal to

- NT i=1 α∈{x,y,z} ∂Ti V ′ i t M α σn αi + σ′ i t M t α V n αi ds. ( 58
)
The free surface boundary condition [START_REF] Benjemaa | Etude de stabilité d'un schéma volumes finis pour les équations de l'élasto-dynamique en maillages non structurés[END_REF] results in NT i=1 k∈K(i) S ik ( V ′ i ) t P ik σ ds = 0. On the other hand, thanks to the absorbing boundary conditions (8), we also have:

1 2 NT i=1 k∈E(i) S ik ( V ′ i ) t P ik σ -A ik V + ( σ′ i ) t P t ik V -B ik σ ds = NT i=1 k∈E(i) S ik ( V ′ i ) t P ik σ + ( σ′ i ) t P t ik V ds, (59)
which closes the proof.

Lemma 4 For any

W h = ( ϕ h , ψ h ) ∈ X 1 h × X 2 h
, the following quantity is positive:

a( W h , W h ) + b( W h , W h ) = 1 2 NT i=1 k∈E(i) S ik ϕ t i A ik ϕ i + ψ t i B ik ψ i ds ≥ 0 . ( 60
) Proof 7 Let be W h = ( ϕ h , ψ h ) ∈ X 1 h × X 2
h a given vector field. Then, we get

a( W h , W h ) + b( W h , W h ) =   NT i=1 α∈{x,y,z} Ti (∂ α ϕ i ) t M α ψ i + ∂ α ψ i t M t α ϕ i dv - 1 2 NT i=1 k∈V(i) S ik ϕ t i P ik ψ i + ψ t i P t ik ϕ i ds - NT i=1 k∈K(i) S ik ψ t i P t ik ϕ i ds - 1 2 NT i=1 k∈E(i) S ik ϕ t i P ik ψ i + ψ t i P t ik ϕ i ds   - 1 2 NT i=1 k∈V(i) S ik ϕ t i P ik ψ k + ψ t i P t ik ϕ k ds + 1 2 NT i=1 k∈E(i) S ik ϕ t i A ik ϕ i + ψ t i B ik ψ i ds .
The first term (between brackets) vanishes by the Green Formula:

Ti (∂ α ϕ i ) t M α ψ i + ∂ α ψ i t M t α ϕ i dv = ∂Ti ϕ t i M α n αi ψ i ds, (61)
combined to the equality ψ t i P t ik ϕ i = ϕ t i P ik ψ i . Moreover, since P ik = α∈{x,y,z} M α n α ik and thanks to the orientation of n ik from T i to T k , we remark that:

NT i=1 k∈V(i) S ik ψ t i P t ik ϕ k ds = - NT i=1 k∈V(i) S ik ϕ t i P ik ψ k ds ,
which implies that the second term is also equal to zero. Therefore, it remains only the third term depending on A ik and B ik , which are respectively a real symmetrical definite positive matrix and a real symmetrical semi-definite positive matrix.

Consequently, if we define the semi-discrete energy by

E h (t) = 1 2 m( W h (t), W h (t)), (62)
where W h (t) = ( V h (t), σh (t)) is the semi-discrete solution of ( 50)-( 51)-( 52), we deduce from Prop. 1 and Lemma 4 that

∂ t E h (t) ≤ 0, ∀t ∈ [0, T ], (63)
which implies that E h (t) is decreasing on [0, T ] and we have

E h (t) ≤ E h (0) ≤ E(0), ∀t ∈ [0, T ]. ( 64
)

Convergence of the semi-discretized problem

We shall at first give a convergence result for the spatially semi-discretized scheme (50)-( 51)-(52). We recall the two following lemmae, well-known in the framework of the finite elements (see [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]):

Lemma 5 Let T i ∈ T h , we assume that u belongs to the space H s+1 (T i ) for s ≥ 0. Let Π be a linear continuous operator from H s+1 (T i ) onto P k (T i ) such that Π(u) = u for all u ∈ P k (T i ). Then we have

u -Π(u) L 2 (Ti) ≤ C 1 h min {s,k}+1 i u H s+1 (Ti) , (65)
u -Π(u) L 2 (∂Ti) ≤ C 1 h min {s,k}+1/2 i u H s+1 (Ti) , (66)
where C 1 is a positive constant only depending on k, s and the regularity parameter ξ of the mesh.

Lemma 6 For all p ∈ P k (T i ), we have

p L 2 (∂Ti) ≤ C 2 h -1/2 i p L 2 (Ti) , (67)
where C 2 is a positive constant only depending on k and the regularity parameter ξ of the mesh.

Thus, we can state the following Lemma:

Lemma 7 Let be W = ( V , σ) the exact solution of ( 6)-( 7)-( 8)-( 9) supposed to belong to C 0 ((0, T ), P H s+1 (Ω) 9 ) and

( V h , σh ) ∈ C 1 ((0, T ), X 1 h × X 2 h
) be the solution of the semi-discrete problem (50)-( 51)-(52). For any t ∈ [0, T ], we have

b( W -P 0 h ( W ), W h -P 0 h ( W ))(t) ≤ Kh min{s,k} ( V h -P 1 h ( V ))(t) 2 L 2 (Ω) 3 + ( σh -P 2 h ( σ))(t) 2 L 2 (Ω) 6 1/2 ( V , σ)(t) P H s+1 (Ω) 9 . ( 68
)
Proof 8 The two previous Lemmae 5 and 6 imply together the following estimates:

S il ( V h -P 1 h ( V )) |T i P il ( σ -P 2 h ( σ)) |T i ds ≤ C P il h min {s,k} i V h -P 1 h ( V ) L 2 (Ti) 3 σ H s+1 (Ti) 6 , (69)
S il ( σh -P 2 h ( σ)) |T i P t il ( V -P 1 h ( V )) |T i ds ≤ C P il h min {s,k} i σh -P 2 h ( σ) L 2 (Ti) 6 V H s+1 (Ti) 3 , (70)
S il

( V h -P 1 h ( V )) |T i A il ( V -P 1 h ( V )) |T i ds ≤ C r + (A il ) h min {s,k} i V h -P 1 h ( V ) L 2 (Ti) 3 V H s+1 (Ti) 3 , (71) S il ( σh -P 2 h ( σ)) |T i B il ( σ -P 2 h ( σ)) |T i ds ≤ C r + (B il ) h min {s,k} i σh -P 2 h ( σ) L 2 (Ti) 6 σ H s+1 (Ti) 6 , ( 72
)
where 71)-(72). In addition, applying (46), we have in the same way, for all l ∈ V(i):

C = C 1 C 2 , l ∈ V(i) ∪ E(i) in (69), l ∈ V(i) ∪ K(i) ∪ E(i) in (70), l ∈ E(i) in (
S il ( V h -P 1 h ( V )) |T i P il ( σ -P 2 h ( σ)) |T l ds ≤ C P il η min{s,k}+1/2 h min {s,k} i V h -P 1 h ( V ) L 2 (Ti) 3 σ H s+1 (T l ) 6 , (73)
S il ( σh -P 2 h ( σ)) |T i P t il ( V -P 1 h ( V )) |T l ds ≤ C P il η min{s,k}+1/2 h min {s,k} i σh -P 2 h ( σ) L 2 (Ti) 6 V H s+1 (T l ) 3 . ( 74
)
Theorem 2 Let (T h ) h be a family of unstructured meshes satisfying (45) and (46). Let ( V , σ) be the exact solution of the symmetrical pseudo-conservative problem ( 6)-( 7)-( 8)-( 9) supposed to belong to C 0 ([0, T ]; P H s+1 (Ω) 9 ) for s ≥ 0, and let

( V h , σh ) ∈ C 1 ([0, T ], X 1 h × X 2 h
) be the solution of the semi-discrete problem (50)-( 51)-(52). Then, there exists a constant K > 0 independent of h such that

max t∈[0,T ] ( V -V h )(t) 2 L 2 (Ω) 3 + ( σ -σh )(t) 2 L 2 (Ω) 6 1/2
≤ KT h min{s,k} ( V , σ) C 0 ((0,T ),P H s+1 (Ω) 9 ) .

(75)

Proof 9 First, applying Lemma 4 to the vector

W h -P 0 h ( W ) = ( V h -P 1 h ( V ), σh -P 2 h ( σ)) ∈ X 1 h × X 2 h , we get a( W h -P 0 h ( W ), W h -P 0 h ( W )) + b( W h -P 0 h ( W ), W h -P 0 h ( W )) ≥ 0, and, then m(∂ t ( W h -P 0 h ( W )), W h -P 0 h ( W )) ≤ m(∂ t ( W h -P 0 h ( W )), W h -P 0 h ( W )) + a( W h -P 0 h ( W ), W h -P 0 h ( W )) + b( W h -P 0 h ( W ), W h -P 0 h ( W )). ( 76
)
Thus, since W h -P 0 h ( W ) ∈ X 1 h × X 2
h , the difference between (56) and (57) gives us

m(∂ t ( W h -W ), W h -P 0 h ( W )) + a( W h -W , W h -P 0 h ( W )) + b( W h -W , W h -P 0 h ( W )) = 0.
After, we substract the previous equality in the right-hand side of (76) and we obtain

m(∂ t ( W h -P 0 h ( W )), W h -P 0 h ( W )) ≤ m(∂ t ( W -P 0 h ( W )), W h -P 0 h ( W )) + a( W -P 0 h ( W ), W h -P 0 h ( W )) + b( W -P 0 h ( W ), W h -P 0 h ( W )). ( 77
)
As P 1 h ( V ) (resp. P 2 h ( σ)) is the orthogonal projection of V (resp. σ) onto X 1 h (resp. X 2 h ) and W h -P 0 h ( W ) ∈ X 1 h ×X 2 h , we necessarily have a( W -P 0 h ( W ), W h -P 0 h ( W )) = 0 (78) and m(∂ t ( W -P 0 h ( W )), W h -P 0 h ( W )) = 1 2 d dt m(( W -P 0 h ( W )), W h -P 0 h ( W )) = 0. ( 79
)
On the other hand, we can check that

m(∂ t ( W h -P 0 h ( W )), W h -P 0 h ( W )) = 1 2 NT i=1 d dt Ti ρ i V h -P 1 h ( V ) 2 (t) + σh -P 2 h ( σ) Λ i σh -P 2 h ( σ) (t) dv (80)
so, by taking also into account (78) and ( 79) in (77), we obtain

d dt ( V h -P 1 h ( V ))(t) 2 L 2 (Ω) 3 + ( σh -P 2 h ( σ))(t) 2 L 2 (Ω) 6 dv ≤ 2 min 1≤i≤NT {ρ i , r -(Λ i )} -1 b( W -P 0 h ( W ), W h -P 0 h ( W ))(t). ( 81
)
Then, we integrate (81) from 0 to t ∈ [0, T ] and, since ( 9) and ( 52) hold, we have

( V h -P 1 h ( V ))(t) 2 L 2 (Ω) 3 + ( σh -P 2 h ( σ))(t) 2 L 2 (Ω) 6 ≤ 2 min 1≤i≤NT {ρ i , r -(Λ i )} -1 T 0 b( W -P 0 h ( W ), W h -P 0 h ( W ))(t)dt. (82)
Therefore, the last inequality becomes

max t∈[0,T ] ( V h -P 1 h ( V ))(t) 2 L 2 (Ω) 3 + ( σh -P 2 h ( σ))(t) 2 L 2 (Ω) 6 (83) ≤ 2 min 1≤i≤NT {ρ i , r -(Λ i )} -1 T 0 b( W -P 0 h ( W ), W h -P 0 h ( W ))(t)dt ( 84
)
which implies, by Lemma 7, that there exists K > 0 such that

max t∈[0,T ] ( V h -P 1 h ( V ))(t) 2 L 2 (Ω) 3 + ( σh -P 2 h ( σ))(t) 2 L 2 (Ω) 6 1/2
≤ KT h min{s,k} ( V , σ) C 0 ((0,T ),P H s+1 (Ω) 9 ) . (85)

On the other hand, we also have by Lemma 5 the existence of C > 0 such that:

∀t ∈ [0, T ], ( V -P 1 h ( V ))(t) 2 L 2 (Ω) 3 + ( σ -P 2 h ( σ))(t) 2 L 2 (Ω) 6 1/2 ≤ Ch min{s,k}+1 V (t) 2 P H s+1 (Ω) 3 + σ(t) 2 P H s+1 (Ω) 6 1/2 , ( 86
)
which combined to (85) concludes the proof.

Convergence of the totally discretized problem

In this section, we prove the convergence of the totally discretized scheme ( 26)-( 27), following the method proposed in [START_REF] Scheid | Convergence of a Discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF].

Theorem 3 Let α be a coefficient such that α = 1 if absorbing boundary conditions are included in the scheme, or α = 2 for domains with only free surface boundary conditions. Let ( V , σ) be the exact solution of ( 6)-( 7)-( 8)-( 9), supposed to belong to

C 2α-1 ([0, T ], L 2 (Ω) 9 ) ∩ C 0 ([0, T ], P H s+1 (Ω) 9 ). Let ( V n h , σn+1/2 h
) be the solution of ( 26)- [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF]. Under a CFL condition of the same type as (38), there exists a constant C > 0 such that the following error estimate holds:

max n=0,...,N V (t n ) -V n h 2 C 0 ([0,T ],L 2 (Ω) 3 ) + σ(t n+1/2 ) - σn+1/2 h 2 C 0 ([0,T ],L 2 (Ω) 6 ) 1/2 ≤ C ∆t α + h min(s,k) ( V , σ) C 2α-1 ([0,T ],L 2 (Ω) 9 ) + ( V , σ) C 0 ([0,T ],P H s+1 (Ω) 9 ) . ( 87
) Proof 10 Let ( V h , σh ) ∈ C 1 ([0, T ], X 1 h × X 2 h
) be the solution of the semi-discrete problem (50)-( 51)-( 52). We are interested in the estimation of the following local consistency error:

ǫ n h = V n h -V h (t n ) 2 L 2 (Ω) 3 + σn+ 1 2 h -σh (t n+ 1 2 ) 2 L 2 (Ω) 6 1/2 , ( 88
)
where

( ϕ h , ψ h ) ∈ X 1 h × X 2 h and ( V n+1 h , σn+1/2 h
) has been computed as follows:

ρ i Ti ϕ t i V n+1 i -V i (t n ) ∆t dv = - α∈{x,y,z} Ti (∂ α ϕ i ) t M α σi (t n+ 1 2 ) dv + k∈V(i) S ik ϕ t i P ik σi (t n+ 1 2 ) + σk (t n+ 1 2 ) 2 ds + k∈E(i) S ik ϕ t i 1 2 P ik σi (t n+ 1 2 ) - 1 2 A ik V i (t n ) ds , ( 89
) Ti ψ t i Λ i σn+ 3 2 i -σi (t n+ 1 2 ) ∆t dv = - α∈{x,y,z} Ti (∂ α ψ i ) t M t α V i (t n+1 )dv + k∈V(i) S ik ψ t i P t ik V i (t n+1 ) + V k (t n+1 ) 2 ds + k∈E(i) S ik ψ t i 1 2 P t ik V i (t n+1 ) - 1 2 B ik σi (t n+ 1 2 ) ds + k∈K(i) S ik ψ t i P t ik V i (t n+1 ) ds . ( 90
)
Using Taylor formulae, there exist c

n+ 1 2 ∈]t n+ 1 2 , t n+1 [ and d n+ 1 2 ∈]t n , t n+ 1 2 [ such that V h (t n+1 ) = V h (t n+ 1 2 ) + ∆t 2 ∂ t V h (t n+ 1 2 ) + ∆t 2 8 ∂ 2 t V h (t n+ 1 2 ) + ∆t 3 48 ∂ 3 t V h (c n+ 1 2 ), (91)
V h (t n ) = V h (t n+ 1 2 ) - ∆t 2 ∂ t V h (t n+ 1 2 ) + ∆t 2 8 ∂ 2 t V h (t n+ 1 2 ) - ∆t 3 48 ∂ 3 t V h (d n+ 1 2 ). ( 92
)
On the other hand, applying (89), ( 91) and (92) in the following equality

Ω ϕ t h V h (t n+1 ) -V n+1 h dv = Ω ϕ t h V h (t n+1 ) -V h (t n ) dv + Ω ϕ t h V h (t n ) -V n+1 h dv, (93)
and, by taking into account (50), there exists e

n+ 1 2 ∈]t n , t n+ 1 2 [ such that Ω ϕ t h V h (t n+1 ) -V n+1 h dv = ∆t 3 48 NT i=1 Ti ϕ t i ∂ 3 t V i (c n+ 1 2 ) + ∂ 3 t V i (d n+ 1 2 ) dv - ∆t 2 4 NT i=1 1 ρ i k∈E(i) S ik ϕ t i A ik ∂ t V i (e n+ 1 2 )ds . ( 94
)
Thus, applying the Cauchy-Schwarz inequality and lemma 6, for domains with only free surface boundary conditions (i.e. E(i) is empty), there exists a constant C > 0 such that:

V h (t n+1 ) -V n+1 h L 2 (Ω) 3 ≤ C∆t 3 V h C 3 ([0,T ],L 2 (Ω) 3 ) , ( 95
)
and when the domain has absorbing boundary conditions, we have the following estimate:

V h (t n+1 ) -V n+1 h L 2 (Ω) 3 ≤ C ∆t 3 V h C 3 ([0,T ],L 2 (Ω) 3 ) + ∆t 2 V h C 1 ([0,T ],L 2 (Ω) 3 ) . ( 96
)
In the very same way, we obtain analogous estimates for σh (t

n+ 3 2 ) - σn+ 3 2 h . Finally, as (∂ 3 t V h , ∂ 3 t σh ) is a discrete approximation of (∂ 3 t V , ∂ 3 t σ)
, the estimation of the local consistency error is bounded in the following way for unbounded domains or domains with free surface boundary conditions (i.e. E(i) is empty):

ε n h ≤ C∆t 3 ( V , σ) C 3 ([0,T ],L 2 (Ω) 9 ) , (97)
and, for domains with absorbing boundary conditions,

ε n h ≤ C ∆t 2 ( V , σ) C 1 ([0,T ],L 2 (Ω) 9 ) + ∆t 3 ( V , σ) C 3 ([0,T ],L 2 (Ω) 9 ) . (98)
For ( ϕ h , ψ h ) ∈ X 1 h × X 2 h , we introduce the notations V n+1 h ( ϕ h ) = NT i=1 V n+1 i ( ϕ i ) = NT i=1 Ti ρ i ϕ t i V i (t n+1 ) -V n+1 i ∆t dv, (99)
S n+ 3 2 h ( ψ h ) = NT i=1 S n+ 3 2 i ( ψ i ) = NT i=1 Ti ψ t i Λ i σi (t n+ 3 2 ) - σn+ 3 2 i ∆t dv, (100)
which, combined with (97) for unbounded domains or domains with free surface boundary conditions, lead to

|V n+1 h | + |S n+ 3 2 h | ≤ C∆t 2 ( V , σ) C 3 ([0,T ],L 2 (Ω) 9 ) ; (101)
using (98) for absorbing boundary conditions, it holds

|V n+1 h | + |S n+ 3 2 h | ≤ C ∆t ( V , σ) C 1 ([0,T ],L 2 (Ω) 9 ) + ∆t 2 ( V , σ) C 3 ([0,T ],L 2 (Ω) 9 ) , (102)
where we denote by |. | the norm of the linear forms on L 2 (Ω) 3 or L 2 (Ω) 6 .

On the other hand, according to (89) and ( 90), the semi-discrete solution

( V h , σh ) ∈ C 1 ([0, T ], X 1 h × X 2 h ) verifies: ρ i Ti ϕ t i V i (t n+1 ) -V i (t n ) ∆t dv = - α∈{x,y,z} Ti (∂ α ϕ i ) t M α σi (t n+ 1 2 ) dv + k∈V(i) S ik ϕ t i P ik σi (t n+ 1 2 ) + σk (t n+ 1 2 ) 2 ds + k∈E(i) S ik ϕ t i 1 2 P ik σi (t n+ 1 2 ) - 1 2 A ik V i (t n ) ds + V n+1 i ( ϕ i ) , ( 103
) Ti ψ t i Λ i σi (t n+ 3 2 ) -σi (t n+ 1 2 ) ∆t dv = - α∈{x,y,z} Ti (∂ α ψ i ) t M t α V i (t n+1 )dv + k∈V(i) S ik ψ t i P t ik V i (t n+1 ) + V k (t n+1 ) 2 ds + k∈E(i) S ik ψ t i 1 2 P t ik V i (t n+1 ) - 1 2 B ik σi (t n+ 1 2 ) ds + k∈K(i) S ik ψ t i P t ik V i (t n+1 ) ds + S n+ 3 2 i ( ψ i ) . (104)
Finally, using ( 26)-( 103) and ( 27)-(104), the field

(V n h , S n+ 1 2 h ) = ( V h (t n ) -V n h , σh (t n+ 1 2 ) - σn+ 1 2 h
) satisfies the 20 following set of equations:

ρ i Ti ϕ t i V n+1 i -V n i ∆t dv = - α∈{x,y,z} Ti (∂ α ϕ i ) t M α S n+ 1 2 i dv + k∈V(i) S ik ϕ t i P ik S n+ 1 2 i + S n+ 1 2 k 2 ds + k∈E(i) S ik ϕ t i 1 2 P ik S n+ 1 2 i - 1 2 A ik V n i ds + V n+1 i ( ϕ i ) , (105)
Ti ψ t i Λ i S n+ 3 2 i -S n+ 1 2 i ∆t dv = - α∈{x,y,z} Ti (∂ α ψ i ) t M t α V n+1 i dv + k∈V(i) S ik ψ t i P t ik V n+1 i + V n+1 k 2 ds + k∈E(i) S ik ψ t i 1 2 P t ik V n+1 i - 1 2 B ik S n+ 1 2 i ds + k∈K(i) S ik ψ t i P t ik V n+1 i ds + S n+ 3 2 i ( ψ i ) . (106)
Now, we introduce the error energy E n at time n∆t by:

E n = 1 2 NT i=1 Ti ρ i (V n i ) t V n i + S n+ 1 2 i t Λ i S n-1 2 i dv - ∆t 8 NT i=1 k∈E(i) S ik (V n i ) t A ik V n i -S n-1 2 i t B ik S n-1 2 i ds. (107)
Reasoning similarly as for the discrete corrected energy, we can prove that under a CFL condition of the same type as (38), we have for domains with free surface boundary conditions only:

V n h 2 L 2 (Ω) 3 + S n+ 1 2 h 2 L 2 (Ω) 6 1/2 ≤ C∆t 2 ( V , σ) C 3 ([0,T ],L 2 (Ω) 9 ) (108)
and, for absorbing boundary conditions, we obtain

V n h 2 L 2 (Ω) 3 + S n+ 1 2 h 2 L 2 (Ω) 6 1/2 ≤ C ∆t ( V , σ) C 1 ([0,T ],L 2 (Ω) 9 ) + ∆t 2 ( V , σ) C 3 ([0,T ],L 2 (Ω) 9 ) . (109)
Finally, we can conclude that the order of convergence of the scheme is O(∆t α + h min(s,k) ) by using the triangular inequality and results already established in the semi-discrete case (see Th. 2):

max n=0,...,N V (t n ) -V n h 2 C 0 ([0,T ],L 2 (Ω) 3 ) + σ(t n+ 1 2 ) - σn+ 1 2 h 2 C 0 ([0,T ],L 2 (Ω) 6 ) 1/2 (110)
≤ C∆t α ( V , σ) C 2α-1 ([0,T ],L 2 (Ω) 9 ) + Ch min(s,k) ( V , σ) C 0 ([0,T ],P H s+1 (Ω) 9 ) (111) ≤ C(∆t α + h min(s,k) ) ( V , σ) C 2α-1 ([0,T ],L 2 (Ω) 9 ) + ( V , σ) C 0 ([0,T ],P H s+1 (Ω) 9 ) , (112)
where α = 1 if the absorbing boundary conditions are invocated, otherwise α = 2 for domains with free surface boundary conditions only.

Numerical results

The numerical method has been applied to two different studies. First, the propagation of an eigenmode in three dimensions of space allows a convergence study of the scheme in the case of uniform and non uniform meshes.

The second application concerns the propagation of the wave produced by an explosive source in a half-space which enables validating source introduction and absorbing boundary conditions and a study of the accuracy of the method in a more realistic context.

Convergence study: propagation of an eigenmode in 3D

The first problem concerns the propagation of an eigenmode in three dimensions of space. The domain of computation is the unit cubic cavity on which we apply free surface boundary conditions. We are interested in the (1, 1, 1) mode whose exact solution at time t and in X = (x, y, z) is given by:

                   V x (t, X) = cos(πx) [sin(πy) -sin(πz)] cos(Ωt) , V y (t, X) = cos(πy) [sin(πz) -sin(πx)] cos(Ωt) , V z (t, X) = cos(πz) [sin(πx) -sin(πy)] cos(Ωt) , σ xx (t, X) = -A sin(πx) [sin(πy) -sin(πz)] sin(Ωt) , σ yy (t, X) = -A sin(πy) [sin(πz) -sin(πx)] sin(Ωt) , σ zz (t, X) = -A sin(πz) [sin(πx) -sin(πy)] sin(Ωt) , σ xy (t, X) = σ xz = σ xz = 0 , (113)
with A = √ 2 ρ µ and Ω = π 2µ/ρ. We set dimensionless values for the medium properties, ρ = 1.0, λ = 0.5 and µ = 0.25, which implies that v p =1. and v s =0.5. The initialisation of the leap-frog scheme is realized by deducing from the analytical expressions (113) the values at t = 0 for the velocity components and at t = ∆t 2 for the stress components. The L 2 -error at step n between the exact value and the solution of the numerical scheme in the unit cube depends on velocities V α (α = x, y, z) at t = n∆t and stresses σ αβ (α, β = x, y, z) at t = (n + 1 2 )∆t and writes

err n L 2 = NT i=1 Ti   α∈{x,y,z} (V α (n∆t, X i ) -(V α ) n i ) 2 +
α,β∈{x,y,z}

σ αβ n + 1 2 ∆t, X i -(σ αβ ) n+1/2 i 2   dv .
Series of different uniform and non uniform meshes have been constructed for a numerical study of the convergence of the scheme. Different methods have been used, the notation P k refering to a scheme based on polynomials of degree k. Uniform meshes are obtained by dividing the domain in cubic cells which are split in six tetrahedra. On the other hand, the non uniform meshes are constructed using the three-dimensional mesh generator GMSH [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with buit-in-pre and post-processing facilities[END_REF] and the different meshes are obtained by successive refinements of an initial coarse mesh. Three of them are presented in figure 1. For both types of meshes, the mesh spacing h is the length of the longest edge of the mesh (i.e. h max ). The characteristics of both uniform and non uniform meshes are given in table 1. On both figures 2 and 3, we display the numerical error, at time t=5.0s, in L 2 -norm and in logarithmic scale, by respect to the mesh spacing h max (left picture) and by respect to the corresponding CPU time (right picture). Figure 2 is concerned with uniform meshes whereas figure 3 refers to non uniform meshes. For uniform and non uniform meshes, we observe a convergence approximatively of order k for a method P k , which proves the enhancement provided by higher degree methods. Note that this study also proves that the free surface condition is accurately discretized, as proved in the previous section. If we examine the efficiency, we notice that, to obtain a given L 2 -error, methods based on high degree polynomials need lower CPU times, for both types of meshes. We also remark the good behavior of the method for non uniform meshes, since CPU times necessary to reach a given error level are comparable to times corresponding to uniform meshes. So, we can conclude that, for this test case, higher order approximations are more accurate and more efficient for both type of meshes.

Explosive source in a half-space

The method has also been applied to a more realistic test case i.e. the propagation of the wave produced by an explosive source in a homogeneous half-space, as illustrated in figure 4, allowing the validation of the source implementation and the absorbing boundary conditions.

The medium is homogeneous of density ρ=2800 kg/m 3 and P-and S-waves velocities are respectively v p =6400 m/s and v s =3700 m/s. We apply a free surface condition on the upper boundary of the domain and absorbing conditions on the other boundaries. The source is located at the center of the domain at 1000m depth (at the coordinates (0;0;1000)) and four surface sensors are placed at 6000m from the epicenter, symmetrically on the xand y-axis, as described in figure 4. The source signal is a Ricker wavelet whose moment-rate is given by

s(t) = 3.10 16 α √ π exp - (t -1) 2 α 2 (114)
with α=0. [START_REF] Moczo | Hybrid modeling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures[END_REF], corresponding to a central frequency equal to f c =1.0Hz and a maximum frequency of f max =3.0Hz. Thus, the wavelength L = v s /f max is approximatively equal to 1233m. The explosive source is introduced as a right hand side on the diagonal components of the stress tensor i.e. for σ xx , σ yy and σ zz . In order to compare the results obtained when using several methods P k , notably the finite volume method (P 0 ), this right hand side writes s(t) g(x, y, z) where g is defined by

g(x, y, z) = 1 M exp - (x -x S ) 2 + (y -y S ) 2 + (z -z S ) 2 h 2 with M = Ω g(x, y, z) dv ,
where (x S , y S , z S ) are the source coordinates. Note that the support of g is taken sufficiently small, compared to the element size, for an accurate approximation of a point source, a right hand side based on a Dirac function s(t) δ(x S , y S , z S ) being only compatible with high degree approximation based methods. Initial conditions for the system are V = 0 and σ = 0 and solutions are calculated until time t=6.0s. The calculation domain is 36km×36km×12km. Solutions have been obtained using several uniform tetrahedral meshes constructed, as previously for the eigenmode test case, by dividing the domain in cubic cells of edge h which are split in six tetrahedra. The mesh spacing h, i.e. here the smallest edge of the mesh, can be expressed as a function of the wavelength L as summarized in table 2 For such a test case where no exact solution is available, we have to define some criteria in order to estimate very meticulously the accuracy of the solution. Firstly, we check two particular properties, consequences of an explosive source, which are i) V x =0 at sensors C1 and C3 (on y-axis) and V y =0 at sensors C2 and C4 (on x-axis) and ii) the evolution at the surface of the vertical velocity component V z only depends on the distance to the epicenter, thus the V z profiles obtained at the four sensors are superimpose. Moreover, for a validation of the source introduction, the amplitudes and the travel times of the computed profiles are compared to a reference solution computed at these four sensors (and not everywhere in the computational domain) using the Discrete Wave Number method [START_REF] Bouchon | Discrete wave-number representation of seismic-source wave fields[END_REF]. Then, three different errors are defined. First, err 1 is concerned with the components which are ideally equal to zero and expresses

err 1 = max t∈[0.;6.] |V x | C1 , |V y | C2 , |V x | C3 , |V y | C4 .
The two other errors depend on the values at the extrema, as illustrated in figure 5. For each non null velocity profile (for instance, here, V z for the reference solution at sensor 1), we identify three extrema and note their values (V z (extr i ), i=1,2,3) and the corresponding times (t i , i=1,2,3).

From these data, we define a relative error, in the amplitude or in time, between two profiles (for instance profiles at sensors C j and C k ), for a given velocity component V α and a given extremum extr i by the following formulae, for amplitudes and times

err rel (V α , C j /C k , extr i ) ampl = V α (extr i ) Cj -V α (extr i ) C k V α (extr i ) ref , err rel (V α , C j /C k , extr i ) t = t i/Cj -t i/C k t i/ref ,
where the subscript ref stands for the corresponding value of the reference solution. Thus, error err 2 measures the asymmetries of the computed solution (i.e. the mean maximum error of the amplitude and the travel time between solutions computed at two sensors) of the vertical velocity component V z which writes

err 2/ampl = 1 3 3 i=1 max 1≤j,k≤4 err rel (V z , C j /C k , extr i ) ampl , err 2/t = 1 3 3 i=1 max 1≤j,k≤4 err rel (V z , C j /C k , extr i ) t .
Finally, err 3 is the mean error between the computed and the reference solutions at all sensors, for all non null velocity components and for all extrema. This error writes, for instance for amplitudes

err 3/ampl = 1 24 3 i=1   j∈{2,4} err rel V x , C j /C j/ref , extr i ampl + j∈{1,3} err rel V y , C j /C j/ref , extr i ampl + 4 j=1 err rel V z , C j /C j/ref , extr i   ,
an equivalent definition being applied to define err 3/t on times. Calculations have been done using methods based on different degrees of interpolation in combination with the meshes defined in table 2 in order to study the convergence and CPU time corresponding to the three types of errors; the results of these calculations are gathered in figure 6 and the CPU time of each calculation is given in table 3. Note that, as previously, P k refers to a method based on polynomials of degree k.

The first line of the figure presents the values of err 1 as a function of the mesh spacing h (left picture) and by respect to the corresponding CPU time (right picture). When examining these results, we notice (left picture) that err 1 corresponding to P 0 method remains nearly constant, whatever the mesh spacing, in the order of 10% of the maximum value of the velocity. Results are clearly improved by the use of higher degree methods since, in comparison with the P 0 results, the value of the P 2 method is reduced by a factor 10 3 . Moreover, higher degree methods P 2 and P 3 are also more efficient since a given error level of accuracy is obtained for lower CPU times. Now, we study the second line of figures corresponding to err 2 , which is a relative error traducing the asymmetry of the V z profiles by respect to h (left picture) and the CPU time (right picture). As previously, we notice that the results of the P 0 method are not accurate enough even using the finest meshes. The improvement on the symmetry is obvious when using the other methods, especially the P 2 and P 3 methods resulting in err 2 values lower than 0.1%. As previously for err 1 , P 2 and P 3 are the two most efficient methods.

Finally, we analyse the results obtained for err 3 (last line of figure 6) which is a relative error between the computed and the reference solutions at the extrema. As previously, the lowest error levels correspond to the P 2 and P 3 methods whereas the P 0 method produces relative errors in the order of 10% and this independently on h. The other methods allow errors lower than 5% when using the finest meshes. Note that, for this last error, if the improvement between P 0 and P 1 is obvious, it is more limited for the two highest degree methods. Furthermore, the P 2 and P 3 methods lead to the same level of error. This fact is also visible when err 3 is plotted as a function of the CPU time (right picture). The improvement on the error levels for the P 2 and P 3 methods is not sufficient to compensate their overcost; thus, for this last error, the P 1 method appears to be most efficient. We also present, in figure 7, the errors err 2/t (left picture) and err 3/t (right picture) on travel times as a function of the mesh spacing. From the study of err 2/t , we firstly notice a clear difference between P 0 and other methods. Relative error on time is approximatively 0.01% for the P 2 and P 3 methods which is negligible, about 0.1% for the P 1 method and greater than 1% for the P 0 method. We remark a comparable behaviour on err 3/t and values about 0.5% for the P 2 and P 3 methods.

In conclusion of this study, we have to distinguish results obtained with the P 0 method from those of the other methods P k (k=1,2,3). First of all, the P 0 method leads to the worst levels of errors, for amplitudes and travel times, and the use of fine meshes do not allow improving these results. Moreover, when combining the different criteria on the solution, the P 2 seems to represent the better compromise between accuracy and efficiency.

The results obtained, at sensor 1, with the method P 2 and a mesh spacing h=200m are presented in figure 8. They are also compared to the reference solution. This constitutes a validation of the source implementation and of the absorbing boundary conditions since no spurious reflections from the boundaries appear in the profiles.

Conclusion

We proposed a discontinuous Galerkin finite element method to solve the first-order hyperbolic system of the elastodynamic equations, written in velocity-stress formulation. This method is applied to non uniform tetraedral meshes allowing an accurate approximation of the medium or mesh refinement in particular areas of the domain. The spatial approximation is based on high degree polynomials defined locally on each element of the mesh. The scheme combines centered fluxes and a leap-frog scheme in time which leads to a non diffusive method. We also detailed a simple absorbing boundary condition, derived from an upwind scheme. In this paper, the properties of the method have been analysed in detail. First, we have proved that this explicit scheme is stable under a CFL type condition resulting from a discrete energy which is preserved in domains with free surfaces and decreasing in presence of absorbing boundary conditions. Moreover, the convergence of the method for the semi-discrete and the fully discrete schemes has been studied. More precisely, for solutions regular enough, the scheme is convergent with an order O(∆t α + h k ), where k is the degree of the polynomial approximation, and α = 2 for domains with only free surface boundary conditions, or α = 1 when the absorbing boundary conditions are included in the scheme. Three-dimensional numerical results have been obtained with a parallel implementation of the DG method. The first one refers to the propagation of an eigenmode in a cubic cavity for which an exact solution is known and a numerical study of the convergence has been realised for series of uniform and non uniform meshes. The second test case concerns a more realistic problem i.e. the propagation of the waves due to an explosion in a half-space involving the validation of the source implementation and the absorbing boundary conditions. The main conclusions of this study are the accuracy and the efficiency of the method particularly for higher degree polynomial interpolations. These results highlight the gain obtained with these methods, especially when compared to the finite volumes method (P0) for which mesh refinement did not enable improving both errors on amplitudes and travel time. The better compromise between accuracy and CPU time seems obtained using the P2 method (based on polynomials of degree 2). Moreover, the free surface condition is accurately approximated. These results are encouraging and could be developed in several interesting directions. First, the treatment at absorbing boundaries could be improved using more accurate boundary conditions. Many techniques have been proposed to reduce spurious reflections, including additional absorbing layers which may be costly in 3D. Then, high-order non reflecting boundary conditions [START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF] or an implicit treatment of our simple absorbing condition could be interesting ways of improvement. Secondly, the global accuracy of the 3D method could be increased by the use of a fourth-order leap-frog scheme in time, as done in [START_REF] Glinsky | A high-order discontinuous Galerkin scheme for elastic wave propagation[END_REF] for 2D elastic wave propagation, leading to high-order methods in both space and time with the ability to use coarser meshes. Finally, we are interested in simulations of more realistic three dimensional problems including topography and realistic media. Thus, the extension of this scheme to viscoelastic rheologies is on the way.

e 2 = e 3 =

 3 -v si , e 4 = e 5 = e 6 = 0, e 7 = e 8 = v si , e 9 = v pi , with

Figure 1 :

 1 Figure 1: Non uniform meshes M2, M3 and M4.

Figure 2 :

 2 Figure 2: Convergence and CPU time for the (1,1,1) eigenmode at t = 5.0s on uniform meshes.

Figure 3 :Figure 4 :

 34 Figure 3: Convergence and CPU time for the (1,1,1) eigenmode at t = 5.0s on non uniform meshes.

Figure 5 :

 5 Figure 5: Velocity V z (m/s) as a function of time (s) ; reference solution at sensor 1. Extrema and corresponding times for the evaluation of errors err 2 and err 3 .

Figure 6 :

 6 Figure 6: Convergence and CPU time for errors err 1 (first line of figures), err 2 (second line of figures) and err 3 (last line of figures) on amplitudes. Explosive source in a half-space, solutions at t =6.s, uniform meshes, computations on 64 processors.

Figure 7 :

 7 Figure7: Convergence of err 2/t (left picture) and err 3/t (right picture) on travel times. Explosive source in a half-space, solutions at t =6.s, uniform meshes, computations on 64 processors.

Figure 8 :

 8 Figure 8: Comparison between computed and reference solutions at sensor 1. V x and V y (left picture) and V z (right picture). Computed solution corresponds to P 2 method for a mesh spacing equal to 200m.

  Therefore, introducing P ik = α∈{x,y,z} M α n α ik , we have:

	α∈{x,y,z} ∂Ti	ϕ t i M α σn αi ds =	k∈V(i) S ik	ϕ t i P ik σds +	k∈E(i) S ik	ϕ t i P ik σds +	k∈K(i) S ik	ϕ t i P ik σds , (14a)
	α∈{x,y,z} ∂Ti

Table 1 :

 1 Characteristics of the uniform and non uniform meshes used for the convergence study ; N N and N T are respectively the number of nodes and tetrahedra of the meshes.

		Uniform meshes				Non uniform meshes
	Mesh	h min	h max	N N	N T	Mesh	h min	h max	N N	N T
	M1	0.5	0.866	27	48	M1	0.707	1.0	14	24
	M2	0.25	0.433	125	384	M2	0.234	0.520	70	227
	M3	0.125	0.216	729	3 072	M3	8.21 10 -2	0.293	398	1 715
	M4	6.25 10 -2	0.108	4 913	24 576	M4	2.90 10 -2	0.193	2 400 11 899
	M5	3.125 10 -2 5.412 10 -2 35 937 196 608	M5	1.03 10 -2 9.52 10 -2 17 993 98 309

Table 2 :

 2 . Calculations have been performed on 64 processors using Correspondance between the mesh spacing h and the wavelength L of the source function. Characteristics of the different meshes (total number of tetrahedra, number of tetrahedra for each processor).Message Passing Interface (MPI), whatever the mesh spacing, for a better comparison of the CPU times between the different methods. The characteristics of the different meshes are also listed in table 2.

	h (in meters)	400	300	200	150	100	75
	Approx. L/n	L/3	L/4	L/6	L/8	L/12	L/16
	Nb. tetrahedra 1 458 000 3 456 000 11 664 000 27 648 000 93 312 000 221 184 000
	Nb. tetra./proc.	22 781	54 000	182 250	432 000	1 458 000	3 456 000

Table 3 :

 3 CPU times of the calculations for the P k methods and the different meshes on 64 processors.

				CPU (s)
	P k /h	400	300	200	150	100	75
	P 0		1min 18s	7min	18min	1h 28min	4h 31min
	P 1		17min	1h 14min	3h 40min 18h 23min 65h 34min
	P 2	38min 1h 44min	8h 12min	25h 20min
	P 3	3h	8h 07min 39h 50 min


Zineb Simeu-Abazi
email: zineb.simeu-abazi@g-scop.grenoble-inp.fr

Maria Di Mascolo
email: maria.di-mascolo@g-scop.grenoble-inp.fr

Eric Gascard
email: eric.gascard@imag.fr

Performance evaluation of centralized maintenance workshop by using Queuing Networks

Keywords: Availability, Maintenance workshop, Performance evaluation, Maintenance cost, Repairing process, Queuing network models

Maintenance workshop integrated in the job-shop system is required in order to maintain the production machines, to ensure the continuity of service and thus to contribute to the improvement of the availability. The objective is to propose a maintenance workshop where all corrective maintenance activities are centralized. When a failed equipment from the production workshop occurs, it is sent to the maintenance workshop in order to be repaired with a given arrival rate. The aim is to maximize the operational availability of the production workshop by reducing the sojourn time in the maintenance workshop. This paper proposes a methodology to design the Central Maintenance Workshop, which enables us to evaluate the performance in terms of cost and sojourn time, for a given budget. For that, we propose a modeling framework based on the queuing network models. Simulation results are given in order to illustrate the influence of different parameters, like arrival rate of the failed equipments and the waiting time of the equipment

INTRODUCTION

The growing complexity of the automation level has influenced almost all the industrial areas. The equipment complexity has forced the companies to increase the productivity performance by optimizing the costs and the manufacturing means. The design of systems that require operational availability cannot be done without a previously developed maintenance structure [START_REF] Monchy | Maintenance : méthodes et organisation[END_REF]. It is the case in some industrial areas, such as aeronautics and the power nuclear stations. The maintenance strategies are applied to the critical systems, whose failures have significant consequences on equipment availability and staff safety. However, in some other areas, the maintenance is not often considered in the system design (Rustenburg, 2000). The cost induced by the failures that cannot be predicted can sometimes damage the entire manufacturing system. Hence, taking into account the maintenance function is essential for the design of a workshop. Obviously, a sophisticated organization and cooperation between the maintenance workshop and the manufacturing system is required. During the development phase, the maintenance concepts are developed. According to the type of the production system, a maintenance workshop is designed and integrated in order to minimize the unavailability times of the failed equipment [START_REF] Daniel | A New Tool: The Generalised and Synchronised Stochastic Petri Nets[END_REF][START_REF] Simeu-Abazi | Maintenance Integration in Manufacturing Systems evaluation and decision International[END_REF]. There are two maintenance types: Preventive maintenance: different actions are programmed in cycle periods; Corrective maintenance: the equipment is repaired after its failure in the maintenance workshop.

Moreover, there is a technique applied in the defense industry, called repair by replacement [START_REF] Fernandez | Framework for implementation of maintenance management in distribution network service providers[END_REF]. The failed equipment is replaced by a new or ready-for-use equivalent one, in order to minimize the waiting time. Concurrently, the failed equipment will be repaired. Obviously, this requires new spare parts in the stock for the initial supply. We propose to integrate this technique in the maintenance workshop. Our aim is to maximize the operational availability of the production workshop by minimizing the sojourn time of equipments in the maintenance workshop. The challenge is to give the structure of the maintenance workshop, performing at a higher level, but under reasonable configuration cost. In order to do this, we develop a method, based on a queuing network model that enables to evaluate the maintenance workshop performance by simulation (Di Mascolo,96). Note that, for relatively simple systems some analytical methods based on the queuing network can be used as in [START_REF] Duri | Performance evaluation and design of base stock systems[END_REF] and [START_REF] Abbou | Influence of the machines failure rate on the maintenance scheduling[END_REF]). When the systems become more complex, simulation must be used. Our aim is to propose a simulation framework for analyzing any queuing system, with an illustration on a maintenance workshop.

The paper is composed of five sections. After the introduction, the second section is dedicated to the presentation of the maintenance workshop and describes its architecture. In the third section, we describe the corrective maintenance methodology. We use the queuing network to model and evaluate by simulation the maintenance workshop performance. Based on the method developed, we determine the repair station to duplicate in order to minimize the equipment sojourn time. In the fourth section, we integrate the technique of repair by replacement. We describe the process, model the functioning and evaluate its performance. The method allows us to determine the stock level.

DESCRIPTION OF MAINTENANCE WORKSHOP

The maintenance workshop includes three principal services [START_REF] Abbou | Influence of the machines failure rate on the maintenance scheduling[END_REF]): 1. Technical service: There is an information exchange between this service and the technical service of the production workshop. The configuration items of the production workshop have been identified in the development phase. The spare parts are selected from the set of configuration items. Then, for each spare part, stock levels have to be determined. Obviously, we assume that the equipment can be transferred into the maintenance workshop. Otherwise, the sectored antennas move to repair the equipment in the production workshop.

The aim is to maximize the operational availability by minimizing the sojourn time of the failed equipment in the maintenance workshop, but under reasonable configuration cost. For that, we develop a method to evaluate the maintenance workshop performances (Simeu-Abazi 1998). We must determine the number of stations included in the maintenance workshop and define the spare parts stock level. From the simulation method developed, the required parameters allow to design suitably the resources of the maintenance workshop, in order to minimize the sojourn time of failed equipment. When the failed equipment arrives in the maintenance workshop, it is put in the input buffer in order to be repaired. The diagnosis determines its failures and describes the repair cycle. The failed equipment is affected to different stations according to failure probabilities. After that, the test is required to validate the repair. Finally, the failed equipment becomes ready-for-use. It is stored in the output buffer before its transfer to the production system. In order to illustrate these aspects, we consider a structure given in figure 1 where different stations are represented: diagnosis, disassembling repair, assembling and test. For each case, we will operate on the adequate parameters.

MODELING AND PERFORMANCE EVALUATION

A method for the evaluation of the maintenance workshop performances is developed by using queuing network formalism. As represented in fig. 1, the maintenance workshop can be composed by a set of interconnected a single service center. The association of one queue to each station enables us to model the maintenance workshop as a network of queues. In each single service center (see model of station in Fig. 4), customers arrive at the service center, wait in the queue if necessary, receive service from the server, and depart. In fact, the service center and its arriving customers constitute a queuing network model with two parameters. First, we must specify the rate at which customers arrive (e.g., one customer every two seconds, or 0.5 customers/second). Second, we must specify the service demand, which is the average service requirement of a customer. These two parameters correspond to the arrival rate of failed equipments denoted by λ and the service rate in each station denoted by µ. For specific parameter values, it is possible to evaluate this model by solving some simple equations, yielding performance measures such as utilization (the proportion of time the server is busy), sojourn time (the average time spent at the service center by a customer, both queuing and receiving service), queue length (the average number of customers at the service center, both waiting and receiving service), and throughput (the rate at which customers pass through the service center) [START_REF] Jackson | Jobshop-like queueing systems[END_REF]. In our application only sojourn time is calculated. It is what characterizes time during which the equipment is unavailable.

We must determine the number of stations included in the maintenance workshop and define the spare parts stock level. From the method developed, the required parameters allow to design suitably the resources of the maintenance workshop, in order to minimize the sojourn time of failed equipment.

Model description

We consider a standard repair process as shown in Fig. 2. The equipment is composed of different parts. We assume that the equipment failure is due to precisely one part, which is detected in the production system. Also, we assume that all failures can be repaired. The entire equipment is transferred The proposed approach is illustrated on an example given in fig. 2. The input of the maintenance workshop corresponds to arrival rate λ of failed equipment. Only two failure modes are considered. Then, according to the diagnostic result obtained by M1, the failed equipment is sent to station M 2 with probability 0.7 for repair 1, or to station M 3 with probability 0.3 for repair 2. Then, the repaired part is sent to adjustment (M4) before assembly in M5. The assembled equipment is sent to the test station M6 for the final verification. For each station, we evaluate the sojourn time w i and the total sojourn time for each equipment denoted by W T . Fig. 2. Example of a repair process.

Modeling with Queuing Network

Queuing networks are used for the modeling of the repair process. For all the equipments, failures occur according to a Poisson process with rate λ (per hour). Each station includes one server with service lead to the time

i i T µ 1 = (hours).
By associating one station to each step of repair process, one obtains the queuing network model shown in figure 3. For a given arrival rate λ, the first block generate parts arrival using intergeneration times from an exponential distribution.

The second block needs some inputs like width arrival, size of queue and the service rate. If the width arrival is upper than one, the traffic merging is generated with a given probabilities. The simulation of this queue single server enables us to obtain the average rate of arrival (to be compared to the input data) as well as the average of waiting time.

Performance evaluation

The performances of the maintenance workshop are evaluated through the sojourn time of the failed parts in this workshop. The objective is to propose a structure of workshop which makes it possible to minimize this time corresponding to the unavailability of the equipment.

From the queuing network model, the maintenance workshop performance can be evaluated.

We assume that the progress of the customers in the stations is Markovian and the service time is exponentially distributed.

For this application, we chose: λ = 0.2/hour, which corresponds to 8 equipments per week and the simulation results are given for an average of 10 000 parts which leads to a relative error of around 0.5% on the service rate.

The obtained results for the configuration of maintenance workshop in figure 3 are given in the following (Table 1).

Improvement of the performance by duplication

We have two methods to reduce the sojourn time. Either we duplicate the stations, which implies an important response lead time, or we store new interchangeable equipments and parts in the stock in order to replace the failed equipments or parts. Therefore, we present here the first method using the same example. The second manner will be explained in the next section.

The simulation result given in table 1, shows that the duration of repair of failed equipment is on average equal to 11.94 h. We notice that the waiting times of equipment in the stations M2 and M5 are the most significant. A solution is to duplicate the stations whose waiting time is important. We duplicate each station M 2 , M 3 and M 5 , one by one and evaluate for each case, the profit for 5 years.

a-Station M 2 duplicated

We duplicate the station M 2 by adding another station. We have the stations M 21 and M 22 , with the same service rate µ 21 = µ 22 =µ 2 . Also, the probability that the failed part will be transferred to the station M 21 or to the station M 22 is 0.35. The queuing network model for the new repair process is shown in Fig. 5. The performances with M 2 duplicated are given in Table 2. The duplication of M2 does not give satisfying results because the average sojourn time decreases only by 0.58 h.

b-Station M 3 duplicated

The queuing network model for the new repair process is shown in Fig. 6. We conserve the same input data. The performances with M 3 duplicated are given in the Table 3.

M1 M2 M 31 M 32 M4 M5 M6

Service rate μ i 0,65 0,4 0,48 0,48 2,00 0,40 1,80

Average sojourn time at station i w i 2.22 3.84 2.22 2.22 0.55 5.00 0.62

Average sojourn time at workshop WT 11.76

The duplication of M3 does not give satisfying results because the average sojourn time decreases only by 0,18 u.t. The duplication of M2 gives the better results.

c-Station M 5 duplicated

The queuing network model for the new repair process is shown in Fig. 7.

By duplicating M5, we have the stations M 51 and M 52 , with the same service rate µ 51 = µ 52 =µ 5 . Also, the probability that the failed part will be transferred to the station M 51 or to the station M 52 is 0,5. The performances with M 5 duplicated are given in Table 4. The duplication of the M5 station makes it possible to minimize the sojourn time in the maintenance workshop.

Compared to the cost of this solution, the investment on a station makes it possible to have a significant profit. Let us suppose that the capital cost of a new station (M2 or M5) is estimated at C Res = 7 500 K €; the cost associated with the unavailability of the equipment is of 500 € /h and the annual operating cost with OC A = 600 €/year. By comparing the results obtained for the 3 cases (the initial model, the model with duplicated M2 and model with duplicated M5) compared to the costs which each case generates, we obtain the results given in table 5. The loss of production and the capital cost are the parameters which make it possible to judge performances of this workshop of maintenance. Then for each structure, the loss of production, the capital cost and the profit at the end of 5 years are calculated to choose the structure of workshop which gives best the performances. The results show that the structure with M5 duplicated gives the best performances.

Table 5 -Performance comparison

Initial model M2 duplicated M5 duplicated

MAINTENANCE OPTIMISATION

The method that enables us to minimize the sojourn time of equipments in the maintenance workshop is the integration of the repair by replacement technique. As it was explained before, interchangeable equipments and parts have to be kept on stock in order to react rapidly to unexpected equipments failures.

The maintenance by replacement can be carried out at various levels. Two levels of replacement are taken into account:

-The replacement of level 1 consists in replacing failed equipment by equipment while drawing from stocks of the repaired equipment in the workshop of maintenance.

-The replacement of level 2 consists to replace part of the failed equipment by another in stocks of the repaired equipment. In this case a station of disassembling is necessary.

We consider now the repair process as described in Fig. 3 with two levels of replacement described above.

Replacement level 1

The principle of this technique is as follows: as soon as faulty equipment is present in the maintenance workshop, its repair starts. In parallel, if another equipment is available in stock, this one is taken for the replacement. Thus, the unavailability time is reduced to its installation and transfer time of the new equipment. The queuing network model corresponding to this workshop with replacement of level 1 is illustrated by the figure 8. For this case, the queuing network model is a base stock system composed of one stage. The purpose is to satisfy requests or demands and to lead the stock of equipments to its nominal level S.

Failures occur according to a Poisson process with rate λ.

When failed equipment arrives at the maintenance workshop, it is immediately sent into repair and, simultaneously, a request for ready-for-use equipment is sent. The failed equipments follow the upper stream of the process described in Fig. 7 and the requests for ready-for-use equipments follow the lower stream.

If equipment is available in the stock of level S, it is picked from the stock to replace the failed equipment in the production workshop, and hence the request is satisfied. In the maintenance workshop, when the equipment is repaired, it becomes a ready-for-use and constitutes an increment for stock S. Table 6 describes different performances of the maintenance workshop according to the stock level of equipments. These results show well that the higher the storage capacity is, the lower sojourn times are. It is of course necessary to make a trade-off between the unavailability cost and the storage cost.

Replacement level 2

The queuing network model of the repair and replacement process is shown in Fig. 9. The queuing network model is a base stock system composed of two stages. In the first stage, we have stock of equipments with nominal level S 2 . In the second one, we have stock of parts with nominal level S 1 .

We assume that: The part price is about 4 K Euro (guarantee 1 year). Annual additional expenses estimated at 2 K Euro. Different combinations from the stock levels are established. Fig. 10 shows, respectively, the investigation costs for the various levels of equipments and parts stocks. That allows carrying out the trade-off between the investigation cost in order to reduce the sojourn time in the maintenance workshop and the benefit on 5 years. As example, for the data used, it is preferable to store 2 parts and 3 equipments ((S1; S2) = (2; 3)), to only store 3 parts ((S1; S2) = (3; 0)). Indeed, starting from the results obtained, for the case (S1; S2) = (2; 3) where the waiting time is equal to 1.74 hour/week, the capital costs during the 5 years are estimated at 124.70 K € and profit at the end of 5 years with 3481.44 K Euros. On the other hand, for the case (S1; S2) = (3; 0), the waiting time is equal to 2.69 lead to the capital costs estimated at 128.70 K€ with a profit at the end of 5 years with 3205.11 K Euros only.

It is obvious that if we investigate in spare parts for the initial supply, we will considerably reduce the waiting time, and hence the operational availability will be increased. However, limited budget is imposed. Then, the stock levels have to be determined.

CONCLUSION

In this paper, we are interested in the maintenance workshop. We have defined this workshop and determine its different services. The corrective maintenance methodology is described. We have proposed to integrate the repair by replacement technique in the maintenance workshop, in order to reduce the waiting time. In order to evaluate the maintenance workshop performances, we have developed a simulation framework based on the Queuing networks. For each situation, the required parameters have allowed reducing the sojourn time, and hence improving the profit. Optimization consists in researching the best trade-off between the availability resulting from the repaired equipment, the size of the maintenance workshop and the maintenance costs. As further directions of our work, we propose to develop an algorithm enabling to determinate automatically the stations that must be duplicated and define the stock level.

Fig. 1

 1 Fig.1Various stations for the repair process.

  workshop. It is dismounted and the failed part is diagnosed.

Fig. 3 .Fig. 4 .

 34 Fig.3. Queuing Network model for repair process.

Fig. 5 .

 5 Fig.5. Queuing Network model with station M2 duplicated.

Fig. 6 .

 6 Fig.6. Queuing Network model with station M3 duplicated.

Fig. 7 .

 7 Fig.7. Queuing Network model with station M5 duplicated.

Fig. 8 .

 8 Fig.8. Queuing Network model for replacement level 1.

Fig. 9 .

 9 Fig.9. Queuing Network model for replacement level 2.Failures occur in the maintenance workshop according to a Poisson process with rate λ (per hour). The failed equipment is sent into repair and, simultaneously, two requests, one for ready-for-use equipment and one for request for ready-foruse part, are sent. If a new or ready-for-use equipment is available in the stock of level S 2 , it replaces the failed equipment. Then, after diagnosis and dismantling in the station M 1 , the defective part is sent into repair. On the same way, when a new or ready-for-use part is available in the stock of level S 1 , it replaces the defective one. The queuing network model is a base stock system composed of two stages. In the first stage, we have stock of equipments with nominal level S 2 . In the second one, we have stock of parts with nominal level S 1 .

Fig. 10

 10 Fig.10 Study of the variation of the sojourn time to the storage capacity of S1 and S2

Table 1 -Simulation values of average sojourn time for each station

 1

	10017 Parts (<λ>=0.2005)	M 1	M 2	M 3	M 4	M 5	M 6
	Service rate	μ	0,65	0.40	0.48	2,00	0,40	1,80
	Average sojourn time at station i	wi	2,22	4,05	2,34	0,54	4,99	0,62
	Average sojourn	W T			11.94
	time at workshop

Table 2 -Simulation results of average sojourn time when station M2 is duplicated

 2

	10043 Parts		M1	M 21	M 22	M3	M4	M5	M6
	(<λ>=0.2010)
	Service rate	μi	0,65	0,4	0,40	0,48	2,00	0,40	1,80
	Average	w i	2.21	3.16	2.87	2.37	0.57	5.11	0.63
	sojourn time at
	Average	WT				11.36
	sojourn time at

Table 3 -Simulation results of average sojourn time when station M3 is duplicated

 3

	10043 Parts
	(<λ>=0.2010)

Table 4 -Simulation results of average sojourn time when station M5 is duplicated.

 4

	10043 Parts		M1	M2	M3	M4	M51	M52	M6
	(<λ>=0.2010)
	Service rate	μi	0,65	0,40	0,48	2,00	0,40	0,40	1,80
	Average sojourn	wi
	time at station i		2.29	3.88	2.35	0.56	3.36	3.37	0.61
	Average sojourn time at workshop	WT			10.28

Table 6 -Simulation results: sojourn time as a function of the size of stock

 6

	Stock size	0	1	2	3	4
	Average sojourn time
	at workshop	11.94 7.25 3.44 1.82 0.93