B Tie

A.-S Mouronval

Systematic development of upwind numerical fluxes for the space discontinuous Galerkin method applied to elastic wave propagation in anisotropic and heterogeneous media with physical interfaces

Keywords: Space discontinuous Galerkin method, Elastic wave propagation, Anisotropy, Heterogeneous medium with physical interfaces, Polycrystalline materials

This research work presents, within a unified and wave oriented variational framework, systematic development of upwind numerical fluxes for the space discontinuous Galerkin methods to model elastic wave propagation in multidimensional anisotropic media with discontinuous material properties. Both first-order velocity-stress and velocity-strain wave formulations are considered. The proposed approach allows the derivation of upwind numerical fluxes in a well structured and hierarchical way according to the degree of inhomogeneity across a physical interface. The numerical fluxes that are exact solutions of a relevant Riemann problem defined at a physical interface are obtained. The developed explicit and intrinsic tensorial expressions of upwind numerical fluxes in multidimensional case allow a better understanding and analysis of the physical meaning of involved terms. As numerical applications, an example with a physical interface separating two materials, one anisotropic and the other isotropic, and an example of polycrystalline material that presents a particular case with a larger number of physical interfaces, are considered. The proposed numerical fluxes are numerically investigated and validated.

Introduction

The space discontinuous Galerkin (dG) method is based on the use of spatially element-wise discontinuous finite element basis functions, and developing appropriate numerical fluxes on element interfaces is a key point for its success [START_REF]Discontinuous Galerkin Methods[END_REF][START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Käser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -I. the two-dimensional isotropic case with external source terms[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][START_REF] De La Puente | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV. anisotropy[END_REF][START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF][START_REF] Ward | A discontinuous Galerkin method for poroelastic wave propagation: The two-dimensional case[END_REF][START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF][START_REF] Zhan | An exact Riemann solver for wave propagation in arbitrary anisotropic elastic media with fluid coupling[END_REF][START_REF] Zhan | Unified Riemann solution for multiphysics coupling: Anisotropic poroelastic/elastic/fluid interfaces[END_REF]. For this purpose, exact solving of the Riemann problem defined on element interfaces is usually recommended. However, when it is applied to elastic media, upwind numerical fluxes solving exactly the Riemann problem can be easily done only in the case of continuous material properties due to the involvement of the fourth order elastic tensor. Recently, research work has been proposed for the derivation of numerical fluxes at physical interfaces, i.e., in the presence of material discontinuities, in 2D [START_REF] Ward | A discontinuous Galerkin method for poroelastic wave propagation: The two-dimensional case[END_REF] or in 3D [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF][START_REF] Zhan | An exact Riemann solver for wave propagation in arbitrary anisotropic elastic media with fluid coupling[END_REF][START_REF] Zhan | Unified Riemann solution for multiphysics coupling: Anisotropic poroelastic/elastic/fluid interfaces[END_REF] cases.

All these works consider only the velocity-strain wave formulation for the elastic media, and they are all placed in the multiphysics context coupling respectively, isotropic elastic/acoustic media for Wilcox et al. [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF], isotropic poroelastic/elastic media for Ward et al. [START_REF] Ward | A discontinuous Galerkin method for poroelastic wave propagation: The two-dimensional case[END_REF], anisotropic elastic media/fluid media for Zhan et al. [START_REF] Zhan | An exact Riemann solver for wave propagation in arbitrary anisotropic elastic media with fluid coupling[END_REF], or anisotropic poroelastic/elastic/fluid interfaces for Zhan et al. [START_REF] Zhan | Unified Riemann solution for multiphysics coupling: Anisotropic poroelastic/elastic/fluid interfaces[END_REF].

Besides, it is noteworthy that penalty fluxes, an alternative to upwind fluxes for dG methods, have also been developed in the literature [START_REF] Warburton | A low-storage curvilinear discontinuous Galerkin method for wave problems[END_REF][START_REF] Ye | A discontinuous Galerkin method with a modified penalty flux for the propagation and scattering of acousto-elastic waves[END_REF][START_REF] Chan | Weight-adjusted discontinuous Galerkin methods: Matrix-valued weights and elastic wave propagation in heterogeneous media[END_REF][START_REF] Shukla | A weight-adjusted discontinuous Galerkin method for the poroelastic wave equation: Penalty fluxes and micro-heterogeneities[END_REF]. To derive numerical fluxes, the penalty method is energy stable and significantly less complicated. The lack of dependence of penalty fluxes on the stiffness matrix allows for a unified and efficient implementation [START_REF] Chan | Weight-adjusted discontinuous Galerkin methods: Matrix-valued weights and elastic wave propagation in heterogeneous media[END_REF]. However, by solving a Riemann problem consistent with physical continuous interface conditions, an upwind flux can be derived, and it is superior, i.e., more accurate, to penalty fluxes. Moreover, such an upwind flux is completely determined by the physical problem itself, and there is no need to deal with the choice of a "numerical" parameter as the penalty parameter for a penalty flux.

The aim of the present work is, by using the unified and wave oriented variational framework that we previously proposed in [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF], to present a systematic development of upwind numerical fluxes in the most general case of multidimensional anisotropic elastic media with discontinuous material properties. Unlike the above-cited works, both first-order velocity-stress formulation defined in [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Käser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -I. the two-dimensional isotropic case with external source terms[END_REF][START_REF] De La Puente | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV. anisotropy[END_REF] and velocity-strain one defined in [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF][START_REF] Ward | A discontinuous Galerkin method for poroelastic wave propagation: The two-dimensional case[END_REF][START_REF] Zhan | An exact Riemann solver for wave propagation in arbitrary anisotropic elastic media with fluid coupling[END_REF] of the elastic wave propagation are considered, and in both cases, closed-form expressions of exact solutions of a Riemann problem defined at a physical interface are obtained. Owing to compact and intrinsic forms expressed in terms of tensors, the derivation of upwind numerical fluxes can be done in a well structured and hierarchical way, which allows better understanding and analysis of the physical meaning of the developed upwind numerical fluxes.

We note that, for the first-order velocity-strain formulation of elastic waves, the numerical flux derived in the present work is the same as the one presented in the above-cited works [START_REF] Ward | A discontinuous Galerkin method for poroelastic wave propagation: The two-dimensional case[END_REF][START_REF] Zhan | An exact Riemann solver for wave propagation in arbitrary anisotropic elastic media with fluid coupling[END_REF]. Our contribution is a new method to derive the upwind flux, using a coordinate-free vector and tensor notation, and exploiting the elastic wave oriented eigenanalysis of a hyperbolic system. The formalism proposed by Wilcox et al. [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF] is close to ours, in the sense that the use of the tensorial forms of strain and stresse fields is preserved, but it is only developed for the isotropic case and for the first-order velocity-strain formulation. While, the formalism used in [START_REF] Ward | A discontinuous Galerkin method for poroelastic wave propagation: The two-dimensional case[END_REF][START_REF] Zhan | An exact Riemann solver for wave propagation in arbitrary anisotropic elastic media with fluid coupling[END_REF][START_REF] Zhan | Unified Riemann solution for multiphysics coupling: Anisotropic poroelastic/elastic/fluid interfaces[END_REF] is a "classical" one and is very different to ours, because the strain and stress tensors are reformulated under a vector form.

As numerical applications, firstly, an example with a physical interface separating two materials, one anisotropic and the other isotropic, is studied. Then, polycrystalline materials that present a particularly interesting case of heterogeneous media with a very large number of physical interfaces are considered. Ultrasonic wave propagation in single-phase and untextured polycristalline materials composed of a large number of elliptic grains are simulated. The proposed numerical fluxes are numerically investigated and validated.

The paper is organized as follows. In Section 2, the unified variational framework using intrinsic tensorial expressions that was proposed in [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF] is recalled for both first-order velocity-stress and velocity-strain hyperbolic systems of elastic wave equations. Section 3 presents the development of different upwind numerical fluxes for the first-order velocity-stress wave system. The necessity of defining a mechanically relevant Riemann problem at a physical interface is particularly discussed. Section 4 applies the proposed approach to the case of the first-order velocity-strain wave system. Finally, in Section 5, the proposed numerical fluxes are applied to the cases where one or more physical interfaces exist.

First-order elastic wave governing equations and their variational frameworks

The governing equations of elastic wave propagation are firstly given in the form of a first-order velocity-stress system with velocity and stress field as primary unknowns [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Käser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -I. the two-dimensional isotropic case with external source terms[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][START_REF] De La Puente | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV. anisotropy[END_REF] (see Eqs. (22.14-22.15) in [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]). Then, the unified and wave oriented variational framework is introduced within which the eigen characteristics of the first-order hyperbolic system are presented and the variational discontinuous Galerkin formulation is given. Finally, the first-order velocity-strain wave equations used by Wilcox et al. [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF] are also considered in the present work, as only isotropic media with physical interfaces were studied by those authors.

First-order velocity-stress elastic wave equations

We consider the wave propagation in an elastic medium Ω ⊂ R d of space dimension d (d = 1, 2, 3) and in a time interval [0, T ]. The first-order velocitystress governing equations can be written as follows: ∀(x, t) ∈ Ω×]0, T [

∂ t U + A ∂x (U ) = 0 or ∂ t v -ρ -1 Div x σ = 0 ∂ t σ -C : ε(v) = 0 (1)
The tensorial compact form in (1) has been proposed by the authors in [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF]. The

generalized unknown U (x, t) = (v(x, t) σ(x, t)) T is composed of v the velocity
unknown and σ the stress unknown, with ( • ) T the adjoint operator. Hence, operator A ∂x and its adjoint A ∂x,T (useful hereafter) are defined as follows:

U (x, t) is a field in R d × R
∀W = (w τ ) T A ∂x w τ = -ρ -1 Div x τ -C : ε(w) , A ∂x,T w τ = -Div x (C : τ ) -ρ -1 ε(w) (2)
where ρ denotes the density, C the fourth-order elasticity tensor, and ":" the usual double dot product between two tensors defined as (C :

ε) ij = C ijkl ε kl .
Herein, the Einstein summation convention is systematically used, and all the vectors and tensors are denoted using bold letters.

The two equations of (1) express respectively the elastodynamic equilibrium and the time derivative of the Hooke's law of elasticity. No body force is considered in the equilibrium equation without loss of generality of the purpose of the present work.

It is useful to recall that, according to the definition of the second-order infinitesimal strain tensor ε, we have:

ε(w) = 1 2 D x w + D T x w = ∂w ∂x i ⊗ s e i (3)
with the following usual space gradient and divergence operators defined using an orthonormal basis (e i ) i=1,...,d :

D x w = ∂w ∂x i ⊗ e i , Div x τ = ∂τ ∂x i • e i (4)
In ( 4), "⊗" denotes the usual tensor product between two vectors: (a ⊗ b) ij = a i b j and " • " the usual dot product between a tensor and a vector:

(A • a) i =
A ij a j . In (3), "⊗ s " is the symmetrized tensor product defined as: (a

⊗ s b) ij = 1 2 (a i b j + a j b i ). Otherwise, it is useful to define the dot product in the vectorial space R d × R d×dsym : ∀W i = (w i τ i ) T , (i = 1, 2), W 1 • W 2 = w 1 • w 2 + τ 1 : τ 2 (5)
According to (4), it is easy to show that, on the boundary ∂D of any subdomain D ⊆ Ω, the flux operator F n (with n = n i e i the outward unit normal vector defined on ∂D) associated to the first-order system (1) is in fact equal to A n , the Jacobian operator in the n direction: ∀W = (w τ ) T ,

F n (W ) = A n (W ) = -ρ -1 τ • n -C : (n ⊗ s w) (6)
In [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF], the subscript index "n" indicates the dependency of F n and of A n on n. In the following, the local orthonormal basis defined on ∂Ω will be denoted by (n, {t α } α=1,...,d-1 ).

Finally, to complete the definition of the elastic wave propagation problem, the following boundary conditions are prescribed:

σ • n = g, on ∂Ω N ×]0, T [ (7a) v = ∂ t u D , on ∂Ω D ×]0, T [ (7b)
The first equation (7a) expresses the Neumann boundary conditions with imposed surface loadings g, and the second one (7b) expresses the Dirichlet boundary conditions with prescribed displacements u D . We recall that the conditions ∂Ω N ∪ ∂Ω D = ∂Ω and ∂Ω N ∩ ∂Ω D = ∅ should be verified. The following initial conditions are also necessary:

σ(x, 0) = C : ε(u 0 (x)), ∀x ∈ Ω (8a) v(x, 0) = v 0 (x), ∀x ∈ Ω (8b)
To develop the variational framework for the space dG method, we recall that an approximated solution U h = (v h σ h ) T of the generalized unknown U = (v σ) T is sought for, and, in contrast to continuous finite element methods, U h is discontinuous from one element to another.

Let M h = {Ω k } k denote a finite element mesh of Ω. In the following, any element Ω k of the mesh M h will be denoted by E and any of its neighboring elements by E . The discontinuous solutions in E and E are respectively denoted by U h and U h . Then, the space dG variational formulation of the elastic wave model problem (1) for any element E can be put into two equivalent forms:

∀W h (x) = (w h (x) τ h (x)) T (W h , ∂ t U h ) E -(A ∂x,T (W h ), U h ) E + < W h , F n (U h , U h ) > ∂E = 0 (9a) (W h , ∂ t U h ) E + (W h , A ∂x (U h )) E + < W h , F n (U h , U h ) -F n (U h ) > ∂E = 0 (9b)
In ( 9), the inner products are defined as follows by recalling the definition of the dot product given in ( 5):

(w 1 , w 2 ) E = E w 1 • w 2 dV, (τ 1 , τ 2 ) E = E τ 1 : τ 2 dV (W 1 , W 2 ) E = E W 1 • W 2 dV = (w 1 , w 2 ) E + (τ 1 , τ 2 ) E (10)
One of the basic ideas of the space dG method is to replace the discontinuous flux F n (U h ) on the element boundary ∂E by a numerical flux F n (U h , U h ).

Hence, an appropriate choice of the numerical flux is essential for the success of the space dG method. Only the definition of the numerical flux on the interior element boundary ∂E int = ∂E\(∂E ∩ ∂Ω), which should also depend on the solution U h in the neighboring elements E of E, is studied in the present work.

As for the numerical flux on external element boundary ∂E ext = ∂E ∩ ∂Ω, a ghost neighbor element E having the same mechanical behavior as E is introduced [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]. In E , it should be imposed that:

σ h • n = 2g -σ h • n and v h = v h for the Neumann boundary conditions on ∂Ω N , and σ h • n = σ h • n and v h = 2∂ t u D -v h for the Dirichlet boundary conditions on ∂Ω D .
Otherwise, we have chosen in our numerical studies to implement the second variational dG formulation (9b), because the terms involved in the differential operator A ∂x are slightly more conventional than those involved in the differential operator A ∂x,T , which is invoked by (9a).

First-order velocity-strain elastic wave equations

For the first-order velocity-strain wave equations used in [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF], the generalized unknown U (x, t) = (v(x, t) ε(x, t)) T is composed of v the velocity unknown and ε the strain unknown. The corresponding tensorial compact form can be written as follows: ∀(x, t) ∈ Ω×]0, T [

M (∂ t U ) + A ∂x (U ) = 0 or ρ∂ t v -Div x (C : ε) = 0 ∂ t ε -C : ε(v) = 0 (11)
In [START_REF] Warburton | A low-storage curvilinear discontinuous Galerkin method for wave problems[END_REF], the infinitesimal strain operator ε( • ) defined by (3) should be distinguished from the primary unknown ε = ε(x, t). The operator M and the space derivative operator A ∂x are defined as follows:

∀W = (w τ ) T M w τ = ρw τ , A ∂x w τ = -Div x (C : τ ) -ε(w) (12)
Then, the flux operator F n associated to the first-order system (11) and the Jacobian operator A n in the n direction verify the following equations:

∀W = (w τ ) T , F n (W ) = A n (W ) = -(C : τ ) • n -n ⊗ s w (13)
The corresponding variational framework reads as:

(W h , M (∂ t U h )) E +(W h , A ∂x (U h )) E + < W h , F n (U h , U h )-F n (U h ) > ∂E = 0 (14)
It is worth noticing that the same notations are used in this paper for both velocity-stress and velocity-strain systems but may have different definitions.

Upwind numerical fluxes for the 1st-order velocity-stress system

In this section, upwind numerical fluxes are firstly developed for the firstorder velocity-stress wave equations. It is emphasized that the unified and elastic wave oriented variational framework previously proposed in [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF] allows a simple, compact and intrinsic expression of the Jacobian operator of the firstorder hyperbolic system (1) in terms of its eigenvalues and eigenmodes. As a consequence, a systematic development of numerical fluxes can be formulated in the multidimensional and general case, i.e. anisotropic and heterogeneous media with physical interfaces. Moreover, it allows a better understanding of the physical meaning of the terms involved in the developed numerical fluxes.

Before giving the definition of the Riemann problem on element interfaces, it is useful to recall some important results given in [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF].

First, the normal Jacobian operator A n (defined by ( 6)) of the first-order velocity-stress wave equations (1) can be decomposed as follows, using its two eigenbases:

A n = λ - n,k R - n,k ⊗ L - n,k + λ + n,k R + n,k ⊗ L + n,k (15)
where

(λ - n,k , R - n,k , L - n,k ) k=qL,{qTα} α=1,d-1 and (λ + n,k , R + n,k , L + n,k ) k=qL,{qTα} α=1,d-1
are respectively the strictly negative and positive eigenvalues and the corre- being respectively the velocity of quasi longitudinal and quasi transverse wave modes propagating in the n direction. The subscript indices "qL" and "qT " respectively refer to terms "quasi longitudinal" and "quasi transverse".

Then, it has been shown that the right and left eigenmodes corresponding to the nonzero eigenvalues of A n are the following:

∀k = qL, {qT α } α=1,••• ,d-1 R ± n,k = w n,k -ρ(z ± n,k ) -1 C : (n ⊗ s w n,k ) , L ± n,k = w n,k -(z ± n,k ) -1 n ⊗ s w n,k (16) with z ± n,k = ρλ ± n,k the acoustic impedance, w n,k = 1 √ 2 γ n,k , (γ n,k ) k=qL,{qTα} α=1,••• ,d-1
the unit eigenvectors of the following usual eigensystem of Γ n the Christoffel tensor:

Γ n • γ n,k = λ 2 n,k γ n,k (17)
The definition of the Christoffel tensor Γ n is recalled in the following:

Γ n • w = ρ -1 (C : (n ⊗ s w)) • n , ∀w (18)
For instance, in the 3D case, there are one quasi longitudinal wave mode and two quasi transverse wave modes. We recall that the word "quasi" means that, in contrast to the isotropic case, we have in the general anisotropic case neither pure longitudinal wave mode verifying γ n,qL n nor pure transverse waves modes verifying γ n,qT ⊥ n.

Furthermore, we note that the right and left eigenmodes corresponding to the zero eigenvalues of A n are not recalled herein, as they are not involved in the decomposition of A n (15) and therefore are not involved in flux terms exchanged on element interfaces.

Among existing works, the numerical flux proposed in [START_REF] Käser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -I. the two-dimensional isotropic case with external source terms[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] uses the same velocity-stress formulation (1), but it is only an approximate solution of the Riemann problem defined on a physical interface across which material properties are discontinuous. However, we have shown in [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF] that, when a numerical flux is not the exact solution of the Riemann problem, it works only when the degree of discontinuity across the physical interface is low. Hereafter, a systematic development of upwind numerical fluxes that are exact solutions of the Riemann problem in the multidimensional case, i.e. 2D and 3D, is proposed.

Otherwise, another important result we have shown in the 1D case is that solving exactly the Riemann problem at element interfaces is not a sufficient condition to get physically sound numerical solutions [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF]. Indeed, different equivalent strong forms of the elastic wave problem give rise to different forms of the Riemann problem and consequently to different interface conditions, which are not all physically coherent. To perform the same analysis in the multidimensional case, two upwind numerical fluxes are developed in the following:

the first one is the solution of the Riemann problem directly defined from the first-order velocity-stress system (1), while the second one is the solution of the Riemann problem modified by taking into account the classical mechanical interface conditions in terms of velocity and stress vector fields.

3.1. Upwind flux directly derived from the first-order v-σ system (1)

We consider the interface of two adjacent elements E and E having respectively (ρ, C, U h ) and (ρ , C , U h ) as densities, elastic moduli and initial states (Figure 1). The Riemann problem defines the states that are results of the propagation of the discontinuity U h -U h . In the following, all the equations are written in the 3D case without loss of generality. In the 3D anisotropic case and by assuming λ - n,qL < λ - n,qT1 < λ - n,qT2 , the Rankine-Hugoniot jump conditions for the first-order velocity-stress system (1) are illustrated in Figure 1 and read as [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF][START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF]:

A n (U h -U a h ) = λ - n,qL (U h -U a h ) (19a) A n (U a h -U b h ) = λ - n,qT1 (U a h -U b h ), A n (U b h -U c h ) = λ - n,qT2 (U b h -U c h ) (19b) A n • U c h + A n • U c h = 0 (19c) A n (U b h -U c h ) = λ - n ,qT2 (U b h -U c h ), A n (U a h -U b h ) = λ - n ,qT1 (U a h -U b h ) (19d) A n (U h -U a h ) = λ - n,qL (U h -U a h ) (19e)
We note that the two outward unit normal vectors of E and E on their interface verify n+n = 0. According to the definition of the eigenvalues and eigenvectors of A n , the following equations hold:

U h -U a h = α qL R - n,qL , U a h -U b h = α qT1 R - n,qT1 , U b h -U c h = α qT2 R - n,qT2 (20a)
U h -U a h = α qL R - n ,qL , U a h -U b h = α qT1 R - n ,qT1 , U b h -U c h = α qT2 R - n ,qT2 (20b)
By respectively adding the three first equations given by (19a)-(19b) and the three last equations given by (19d)-(19e) and by applying [START_REF] Tie | Adaptive time discontinuous Galerkin method for numerical modeling of wave propagation in shell and 3D structures[END_REF], the Rankine-Hugoniot jump conditions finally become:

A n (U h ) -A n (U c h ) = α k λ - n,k R - n,k (21a)
A n (U c h ) + A n (U c h ) = 0 (21b) A n (U h ) -A n (U c h ) = α k λ - n ,k R - n ,k (21c)
In the case where

λ - n,qT1 = λ - n,qT2 (resp. λ - n ,qT1 = λ - n ,qT2
), the two arrows λ - n,qT1 and λ - n,qT2 (resp. λ - n ,qT1 and λ - n ,qT2 ) in Figure 1 overlap, and the state

U b h (resp. U b h ) is removed.
Then, the equations (19b) and (20a) become:

A n (U a h -U c h ) = λ - n,qT1 (U a h -U c h ) (22a) U h -U a h = α qL R - n,qL , U a h -U c h = α qT1 R - n,qT1 + α qT2 R - n,qT2 (22b)
which finally lead to the same Rankine-Hugoniot jump conditions (21a). Same remark holds for the quantities related to the element E and the corresponding equations (19d), (20b) and (21c).

Solving the Riemann problem (21) leads to the determination of the six unknown states {U a , U b , U c , U a , U b , U c }, i.e., the six characteristic coefficients {α k , α k } k=qL,qT1,qT2 . Then, upwind numerical fluxes defined as

F n (U h , U h ) = A n (U c h ) = 1 2 (A n (U c h )-A n (U c h )
) can be calculated by using the following two equivalent forms:

F n (U h , U h ) = 1 2 (A n (U h ) -A n (U h ) -α k λ - n,k R - n,k + α k λ - n,k R - n,k ) (23a) = A n (U h ) -α k λ - n,k R - n,k (23b)
It is worth noticing that the two forms (23a) and (23b) are equivalent only in the case where the Riemann problem ( 21) is exactly solved. Otherwise, according to (21b), the flux defined by ( 23) is conservative, i.e., from the point of view of the neighboring element E , we have

F n (U h , U h ) = -F n (U h , U h ).
Now, the solving of the Riemann problem is considered. Two operators ℘ vect ( • ) and ℘ tens ( • ) are defined: when applied to a generalized vector W = (w τ ) T , they respectively give its vectorial and tensorial components:

℘ vect (W ) = w , ℘ tens (W ) = τ (24)
Using ℘ vect and ℘ tens , the vectorial part and the tensorial parts of the equations (21b) can be separately manipulated. According to the definition of the Jacobian operator ( 6) and the definition (18) of the Christoffel tensor Γ n , (21b) is rewritten as follows:

σ c h • n ρ + σ c h • n ρ = 0 , ρΓ n • v c h -ρ Γ n • v c h = 0 (25)
The second equation of ( 25) is obtained by applying the operator n • ℘ tens to Eq.(21b).

Then, by eliminating the two unknown states U c h and U c h in (21) and using [START_REF] Stanke | A unified theory for elastic wave-propagation in polycrystalline materials[END_REF], the following result is obtained.

Theorem 3.1. The characteristic coefficients {α k , α k } k=qL,qT1,qT2 of the Riemann problem [START_REF] Leclère | Parallel and adaptive FE modeling of elastic wave propagation in structures (Modélisation parallèle de la propagation d'ondes dans les structures par éléments finis adaptatifs)[END_REF] are the solution of the following system of linear equations:

[I d ] [B] [B ] [I d ] • {α k } {α k } = { * L - n,k • U h - * * L - n,k • U h } { * L - n ,k • U h - * * L - n ,k • U h } (26)
In

B kl,k =l = C - z,k 2 ρ ρ z - n ,l δz - n ,kl z - n,k z - n ,k γ n,k • γ n ,l , B kl,k =l = C - z,k 2 ρ ρ z - n,l δz - n,kl z - n ,k z - n,k γ n ,k • γ n,l (27)
with δz - n,kl = z - n,k -z - n,l , δz - n ,kl = z - n ,k -z - n ,l
, and :

C - z,k = z - n,k R z - n,k = z - n,k z - n,k V > 0 , C - z,k = z - n,k R z - n ,k = z - n,k z - n,k V > 0 (28) z - n,k R and z - n,k
V respectively denote the harmonic and arithmetic means of the acoustic impedances of the k-th eigenvector.

{ * L - n,k , * * L - n,k , * L - n ,k , * * L - n ,k } are
the perturbed left eigenmodes of {A n , A n } calculated by using the material properties of the adjacent element in the following way:

* L - n,k =   C - z,k ℘ vect (L - n,k ) C - z,k ℘ tens (L - n,k )   , * * L - n,k =      ρρ (z - n,k ) 2 Γ n • ℘ vect ( * L - n,k ) ρ ρ ℘ tens ( * L - n,k )      * L - n ,k =    C - z,k ℘ vect (L - n ,k ) C - z,k ℘ tens (L - n ,k )    , * * L - n ,k =       ρ ρ (z - n ,k ) 2 Γ n • ℘ vect ( * L - n ,k ) ρ ρ ℘ tens ( * L - n ,k )       (29)
Proof. See Appendix A1.

It is worth noticing that the perturbed left eigenmodes

{ * L - n,k , * * L - n,k } and { * L - n ,k , * * L - n ,k
} take into account the coupling between the wave modes of the same type, e.g., both qL-modes, of the adjacent elements E and E . In contrast, the matrices [B] and [B ] take into account the coupling between two wave modes of different types, e.g., a qL wave mode from E with a qT wave mode from E .

We also remark that ( 27) and ( 29) are completely symmetric with respect to E and E .

Before giving one of the main results of the present work, it is worth recalling the numerical flux proposed by Käser et al. [START_REF] Käser | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -I. the two-dimensional isotropic case with external source terms[END_REF], which is expressed only in terms of material properties from the interior of the element E from one side of element interface, i.e.,

F Käser n (U h , U h ) = λ + n,k (R + n,k ⊗ L + n,k ) • U h + λ - n,k (R - n,k ⊗ L - n,k ) • U h (30)
Indeed, it would be interesting to present new terms brought by the upwind numerical fluxes proposed in the present work.

To obtain the upwind numerical fluxes, the linear system [START_REF] Bai | Finite element modeling of ultrasonic wave propagation in polycrystalline materials[END_REF] needs to be solved. The 2d × 2d matrix of the system of ( 26

F n (U h , U h ) = F Käser n (U h , U h ) + λ - n,k (R - n,k ⊗ (L - n,k - * L - n,k )) • U h -(R - n,k ⊗ (L - n,k - * * L - n,k )) • U h (31) 2) 2nd-order upwind numerical flux by taking [R] -1 ≈ [I 2d ] -[∆R] and by remarking that [R].([I 2d ] -[∆R]) -[I 2d ] = [∆R] 2 F n (U h , U h ) = F 1st-order n (U h , U h ) -λ - n,k l =k B kl (R - n,k ⊗ * * L - n ,l ) • U h -(R - n,k ⊗ * L - n ,l ) • U h (32)
where F 1st-order n (U h , U h ) denotes the flux defined by (31).

3) Upwind numerical flux being the exact solution of the Riemann problem(21)

F n (U h , U h ) = F 1st-order n (U h , U h ) -λ - n,k D kl (R - n,k ⊗ * L - n,l ) • U h -(R - n,k ⊗ * * L - n,l ) • U h -λ - n,k H kl (R - n,k ⊗ * * L - n ,l ) • U h -(R - n,k ⊗ * L - n ,l ) • U h (
[R] -1 = [I 2d ] + [D] -[H] -[H ] [D ] (34)
Proof. It is straightforward by using the decomposition (15) of A n and by substituting the solution of {α k } in the asymmetric form of numerical flux (23b).

For several specific case, the following corollaries are immediate from Theorem 3.2:

Corollary 3.2.1. On an element interface with continuous material properties,

we have

[B] = [B ] = 0, * L - n,k = * * L - n,k = L - n,k and * L - n ,k = * * L - n ,k = L - n ,k , so
the characteristic coefficients {α k , α k } k=qL,qT1,qT2 of the Riemann problem [START_REF] Leclère | Parallel and adaptive FE modeling of elastic wave propagation in structures (Modélisation parallèle de la propagation d'ondes dans les structures par éléments finis adaptatifs)[END_REF] can be explicitly solved:

α k = L - n,k • (U h -U h ) , α k = L - n ,k • (U h -U h ) (35)
All the three fluxes defined in Theorem 3.2 are identical. They are all equal to the Käser's flux [START_REF] Tie | Some comparisons and analyses of time or space discontinuous Galerkin methods applied to elastic wave propagation in anisotropic and heterogeneous media[END_REF] and can be written as follows:

F n (U h , U h ) = F Käser n (U h , U h ) (36a) = λ + n,k (R + n,k ⊗ L + n,k ) • U h -λ + n ,k (R + n ,k ⊗ L + n ,k ) • U h (36b)
Proof. (36a) is straightforward. The symmetric form of the numerical flux (36b) can be obtained by remarking that, in the case of continuous material properties, the following equations hold (see Appendix A2):

λ n ,k = λ n,k , R ∓ n ,k = -R ± n,k , L ∓ n ,k = -L ± n,k (37)
Corollary 3.2.2. On an element interface with discontinuous but isotropic material properties on its two sides, we have [B] = [B ] = 0, so the characteristic coefficients {α k , α k } k=qL,qT1,qT2 of the Riemann problem ( 21) can be explicitly solved:

α k = * L - n,k • U h - * * L - n,k • U h , α k = * L - n ,k • U h - * * L - n ,k • U h (38)
All the three fluxes defined in Theorem 3.2 are identical, and the 1st-order flux

F 1st-order n (U h , U h ) defined by (31)
is the exact solution of the Riemann problem [START_REF] Leclère | Parallel and adaptive FE modeling of elastic wave propagation in structures (Modélisation parallèle de la propagation d'ondes dans les structures par éléments finis adaptatifs)[END_REF].

Proof. It is straightforward by simply recalling that in the isotropic case γ n,qL = n, γ n,qT1 = t 1 and γ n,qT2 = t 2

In the case of media with discontinuous anisotropic material properties, the sub-matrices [B] and [B ] of [∆R] defined by ( 27) are proportional to the degree of discontinuity across a physical interface. Thus, when the degree of discontinuity is sufficiently low, the 2nd-order flux, even the 1st-order flux, would be a sufficiently good numerical flux. However, when the degree of discontinuity is high, it is necessary to use the numerical flux that is the exact solution of the Riemann problem. It is worth noticing that the Käser's numerical flux [START_REF] Tie | Some comparisons and analyses of time or space discontinuous Galerkin methods applied to elastic wave propagation in anisotropic and heterogeneous media[END_REF] should be worse than the 1st-order flux (31), as it does not use the perturbed

left eigenmodes { * L - n,k , * * L - n,k }.
Finally, the most important remark concerning the numerical fluxes defined in Theorem 3.2 is that they are not physically sound in the case of media with discontinuous material properties. Indeed, on a physical interface, i.e. ρ = ρ or/and Γ n = Γ n , the interface conditions (25) derived from the equation (19c)

of the Rankine-Hugoniot jump conditions [START_REF] Li | Implementation and adaptivity of space-time finite element method for structural dynamics[END_REF] are not equivalent to the classical interface conditions of continuous velocity and stress vector. We recall that it has been shown in the 1D case [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF] that the so-defined numerical flux was not able to give rise to physically sound wave propagation solutions.

Furthermore, the following result can be proved for the numerical flux (33) that is the exact solution of the Riemann problem:

Theorem 3.3. The numerical flux (33) directly derived from the first-order velocity-stress system (1) is not consistent.

Proof. To prove that a numerical flux F n (U h , U h ) is consistent, it should to be shown that:

F n (U , U ) = F n (U ) (39)
According to (23b), it is necessary to have, when U h = U h = U and for k = qT, qT 1 , qT 2 , α n,k = 0. But, this does not hold in the general case, as the second member of the system of linear equations [START_REF] Bai | Finite element modeling of ultrasonic wave propagation in polycrystalline materials[END_REF], becoming in this

case ({( * L - n,k - * * L - n,k ) • U }, {( * L - n ,k - * * L - n ,k ) • U })
, generally does not equal to zero. In other words, if α n,k = 0 and α n,k = 0 for k = qT, qT 1 , qT 2 with [START_REF] Leclère | Parallel and adaptive FE modeling of elastic wave propagation in structures (Modélisation parallèle de la propagation d'ondes dans les structures par éléments finis adaptatifs)[END_REF] then results in the following equation:

U h = U h = U ,
A n (U | E ) + A n (U | E ) = 0 (40)
that is equivalent to:

σ| E • n ρ + σ| E • n ρ = 0 , ρΓ n • v| E -ρ Γ n • v| E = 0 (41)
But (41) is generally not true when the element interface is also a physical interface.

Upwind flux based on the mechanical interface conditions

In order to remove the physical inconsistency of the interface conditions [START_REF] Stanke | A unified theory for elastic wave-propagation in polycrystalline materials[END_REF], they are modified in the following way to take into account the classical interface conditions, i.e. velocity and stress vector continuities:

σ c h • n + σ c h • n = 0 , v c h -v c h = 0 (42)
Both interface conditions ( 25) and ( 42) are equivalent only in the case where the elastic moduli and the density of the propagating medium are continuous across the interface.

It is easy to show that the interface conditions (42) is in fact a consequence of the following Rankine-Hugoniot jump conditions:

Ãn (U h ) -Ãn (U c h ) = αk λ - n,k M (R - n,k ) (43a) Ãn (U c h ) + Ã n (U c h ) = 0 (43b) Ã n (U h ) -Ã n (U c h ) = α k λ - n ,k M (R - n ,k ) (43c)
which correspond to another equivalent form of the first-order velocity-stress system (1):

M (∂ t U ) + Ã∂x (U ) = 0 or ρ∂ t v -Div x σ = 0 C -1 : ∂ t σ -ε(v) = 0 (44)
with the following definitions of operators:

M w τ = ρw C -1 : τ , Ã∂x w τ = -Div x τ -ε(w) , Ãn w τ = -τ • n -n ⊗ s w (45)
We note that the following equations hold:

M (R n,k ) = ρL n,k , Ãn = M • A n = z ± n,k L ± n,k ⊗ L ± n,k (46)
When the characteristic coefficients {α k , α k } k=qL,qT1,qT2 are calculated, the upwind numerical fluxes to be used for the first-order velocity-stress strong form (44) can be defined in the same way as in the preceding section, by taking into account (43a):

F n (U h , U h ) = Ãn (U c h ) = Ãn (U h ) -αk z - n,k L - n,k (47)
However, it is also possible to keep using the variational formulation [START_REF] Zhan | An exact Riemann solver for wave propagation in arbitrary anisotropic elastic media with fluid coupling[END_REF] of the first-order velocity-stress strong form (1) instead of the corresponding variational formulation of (44). In this case, by taking into account (46), the equation to define upwind numerical fluxes reads as:

F n (U h , U h ) = M -1 • F n (U h , U h ) = A n (U h ) -αk λ - n,k R - n,k (48)
It should be noted that both choices, the strong form (1) with the numerical flux (47) or the strong form (44) with the numerical flux (48), are equivalent, and their numerical implementations give rise to identical results, which is confirmed by our numerical investigation. Moreover, we notice that

F n (U h , U h ) = -F n (U h , U h ), even if the numerical flux (47) is conservative, i.e., F n (U h , U h ) = -F n (U h , U h ).
To solve the Riemann problem (43), the same technique presented in the preceding section is used, and the following result concerning the mechanically based Riemann problem can be proved:

Theorem 3.4. The characteristic coefficients {α k , α k } k=qL,qT1,qT2 of the Riemann problem (43) are the solution of the following linear system of equations:

[I d ] [ B] [ B ] [I d ] • {α k } {α k } = {˜L - n,k • (U h -U h )} {˜L - n ,k • (U h -U h )} (
Bkl,k =l = - C - z,k 2 δz - n ,kl z - n ,k γ n,k • γ n ,l , B kl,k =l = - C - z,k 2 δz - n,kl z - n,k γ n ,k • γ n,l (50)
where

C - z,k , C - z,k , δz - n ,kl and δz - n,kl are defined in Theorem 3.1. {˜L - n,k ,˜L - n ,k }
are the perturbed left eigenmodes of {A n , A n } calculated by using the material properties of the adjacent element in the following way:

˜L- n,k = C - z,k ℘ vect (L - n,k ) C - z,k ℘ tens (L - n,k ) , ˜L- n ,k =   C - z,k ℘ vect (L - n ,k ) C - z,k ℘ tens (L - n ,k )   (51)
Proof. See Appendix A3.

As in the preceding section, by decomposing the 2d×2d matrix of the system of linear equations ( 49

[ R].[I 2d ] -[I 2d ] = [∆ R] F n (U h , U h ) = M ( F Käser n (U h , U h )) + z - n,k (L - n,k ⊗ (L - n,k -˜L - n,k )) • (U h -U h ) (52a) F n (U h , U h ) = F Käser n (U h , U h ) + λ - n,k (R - n,k ⊗ (L - n,k -˜L - n,k )) • (U h -U h ) (52b)
2) 2nd-order upwind numerical flux by taking

[ R] -1 ≈ [I 2d ] -[∆ R] and by remarking that [ R].([I 2d ] -[∆ R]) -[I 2d ] = [∆ R] 2 F n (U h , U h ) = F 1st-order n (U h , U h ) -z - n,k l =k Bkl (L - n,k ⊗˜L - n ,l ) • (U h -U h ) (53a) F n (U h , U h ) = F 1st-order n (U h , U h ) -λ - n,k l =k Bkl (R - n,k ⊗˜L - n ,l ) • (U h -U h ) (53b)
where F 1st-order n (U h , U h ) and F 1st-order n (U h , U h ) denote the fluxes defined respectively by (52a) and (52b).

3) Upwind numerical flux being the exact solution of the Riemann problem(43) As in the preceding section, the following corollaries can be proved for several specific cases:

F n (U h , U h ) = F 1st-order n (U h , U h ) -z - n,k L - n,k ⊗ Dkl ˜L- n,l + Hkl ˜L- n ,l • (U h -U h ) (54a) F n (U h , U h ) = F 1st-order n (U h , U h ) -λ - n,k R - n,k ⊗ Dkl ˜L- n,l + Hkl ˜L- n ,l • (U h -U h ) (54b
αk = L - n,k • (U h -U h ) , α k = L - n ,k • (U h -U h ) (56)
All the three fluxes defined in Theorem 3.5 are identical. They are equivalent to the one defined by (36) in Corollary 3.2.1.

Corollary 3.5.2. On an element interface with discontinuous but isotropic material properties on its two sides, we have [ B] = [ B ] = 0, so the characteristic coefficients {α k , α k } k=qL,qT1,qT2 of the Riemann problem (43) can be explicitly solved:

αk = ˜L- n,k • (U h -U h ) , α k = ˜L- n ,k • (U h -U h ) (57)
All the three fluxes defined in Theorem 3.5 are identical, and the 1st-order flux

F 1st-order n (U h , U h ) defined by (52a)
is the exact solution of the Riemann problem (43).

Concerning the consistency of the numerical flux (54) that is the exact solution of the Riemann problem, the following result can be proved.

Theorem 3.6. The numerical flux (54) based on the mechanical interface conditions and corresponding to the first-order velocity-stress system (44) is consistent.

Proof. As the second member of the system of linear equations ( 49) is equal to zero with U h = U h = U , hence, we have αn,k = α n,k = 0 for k = qT, qT 1 , qT 2 .

Then, according to (47) and (48), we get F n (U , U ) = F n (U ) and F n (U , U ) =

F n (U ).

Upwind numerical fluxes for the first-order velocity-strain system

In this section, we apply the approach developed in the preceding section to the first-order velocity-strain wave equations defined in Section 2.2. In [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF],

an upwind numerical flux that is exact solution of the Riemann problem was developed in the case of isotropic heterogeneous media with discontinuous elastic moduli and density. Hereafter, the more general case of anisotropic media is considered. Firstly, the elastic wave oriented eigenanalysis is applied to the firstorder velocity-strain waves equations [START_REF] Warburton | A low-storage curvilinear discontinuous Galerkin method for wave problems[END_REF]. Secondly, the Riemann problem is defined and solved to develop upwind numerical fluxes.

4.1. Elastic wave oriented eigenanalysis of the first-order v-ε system [START_REF] Warburton | A low-storage curvilinear discontinuous Galerkin method for wave problems[END_REF] As in the case of the first-order velocity-stress system, before giving the definition of the Riemann problem at an element interface, eigenanalysis of the first-order velocity-strain system (11) is performed, i.e., the solving of the following eigenproblem:

A n (R n,k ) = λ n,k M (R n,k ) , A T n (L n,k ) = λ n,k M (L n,k ) (58)
with A n and M defined respectively in [START_REF] Chan | Weight-adjusted discontinuous Galerkin methods: Matrix-valued weights and elastic wave propagation in heterogeneous media[END_REF] and [START_REF] Ye | A discontinuous Galerkin method with a modified penalty flux for the propagation and scattering of acousto-elastic waves[END_REF]. The following results can be proved.

Among

= qL, {qT α } α=1,••• ,d-1 R ± n,k = w n,k -ρ(z ± n,k ) -1 n ⊗ s w n,k , L ± n,k = ρ -1 w n,k -ρ(z ± n,k ) -1 C : (n ⊗ s w n,k ) (59
) As in the previous case of the first-order velocity-stress system,

w n,k = 1 √ 2 γ n,k , with (γ n,k ) k=qL,{qTα} α=1,••• ,d-1 the unit eigenvectors of the Christoffel tensor Γ n .
It can be shown that the following orthonormality relation holds between the two right and left eigenbases:

R n,k • M (L n,l ) = δ kl , ∀k, l = 1, • • • , m (60)
where δ kl denotes the usual Kronecker delta. Hence, using its two eigenbases, the normal Jacobian operator A n defined by ( 13) of the first-order velocity-strain wave equations (11) can be decomposed as follows:

A n = λ - n,k M (R - n,k ) ⊗ M (L - n,k ) + λ + n,k M (R + n,k ) ⊗ M (L + n,k ) (61)

Upwind numerical fluxes

In the 3D case, the Riemman problem corresponding to the first-order velocitystrain equations [START_REF] Warburton | A low-storage curvilinear discontinuous Galerkin method for wave problems[END_REF], i.e., the Rankine-Hugoniot jump conditions across an element interface (Figure 1), reads as:

A n (U h ) -A n (U c h ) = α k λ - n,k M (R - n,k ) (62a) A n (U c h ) + A n (U c h ) = 0 (62b) A n (U h ) -A n (U c h ) = α k λ - n ,k M (R - n ,k ) (62c)
(62b) gives rise to the following interface conditions that are clearly physically sound:

(C : ε c h ) • n + (C : ε c h ) • n = 0 , v c h -v c h = 0 (63)
When the characteristic coefficients {α k , α k } k=qL,qT1,qT2 are calculated, upwind numerical fluxes can be defined in the same way as in the preceding section, by taking into account (62a):

F n (U h , U h ) = A n (U c h ) = A n (U h ) -α n,k λ - n,k M (R - n,k ) (64)
When the Riemann problem is exactly solved by the numerical flux, the latter is conservative, as we have, due to (62b

), F n (U h , U h ) = -F n (U h , U h ).
To solve the Riemann problem (62), the same technique presented in the preceding section is used, and the following result can be proved:

Theorem 4.1. The characteristic coefficients {α k , α k } k=qL,qT1,qT2 of the Riemann problem (62) are the solution of the following linear system of equations:

[I d ] [ B] [ B ] [I d ] • {α k } {α k } = {M ( * L - n,k ) • (U h -U h )} {M ( * L - n ,k ) • (U h -U h )} (
L - n,k = C - z,k ℘ vect (L - n,k ) C - z,k ℘ tens (L - n,k ) , * L - n ,k =   C - z,k ℘ vect (L - n ,k ) C - z,k ℘ tens (L - n ,k )   (66)
Proof. See Appendix A4.

We remark that, on the one side, both systems of linear equations (49) of

[ R].[I 2d ] -[I 2d ] = [∆ R] F n (U h , U h ) = λ + n,k (M (R + n,k ) ⊗ M (L + n,k )) • U h + λ - n,k (M (R - n,k ) ⊗ M (L - n,k )) • U h + λ - n,k (M (R - n,k ) ⊗ M (L - n,k - * L - n,k )) • (U h -U h ) (67) 2) 2nd-order upwind numerical flux by taking [ R] -1 ≈ [I 2d ] -[∆ R] and by remarking that [ R].([I 2d ] -[∆ R]) -[I 2d ] = [∆ R] 2 F n (U h , U h ) = F 1st-order n (U h , U h )-λ - n,k l =k Bkl (M (R - n,k )⊗M ( * L - n ,l )) • (U h -U h ) ( 68
)
where F 1st-order n (U h , U h ) denotes the flux defined by (67).

3) Upwind numerical flux being the exact solution of the Riemann problem(62) As in the preceding section, the following corollaries can be proved for several specific cases: Corollary 4.2.1. On an element interface with continuous material properties,

F n (U h , U h ) = F 1st-order n (U h , U h ) -λ - n,k M (R - n,k ) ⊗ Dkl M ( * L - n,l ) + Hkl M ( * L - n ,l ) • (U h -U h ) (69) with [ D], [ D ], [ H]
we have [ B] = [ B ] = 0, * L - n,k = L - n,k and * L - n ,k = L - n ,k
, so the characteristic coefficients {α k , α k } k=qL,qT1,qT2 of the Riemann problem (62) can be explicitly solved:

α k = M (L - n,k ) • (U h -U h ) , α k = M (L - n ,k ) • (U h -U h ) ( 70
)
All the three fluxes defined in Theorem 4.2 are identical and can be written as follows:

F n (U h , U h ) = λ + n,k (M (R + n,k ) ⊗ M (L + n,k )) • U h + λ - n,k (M (R - n,k ) ⊗ M (L - n,k )) • U h (71a) = λ + n,k (M (R + n,k ) ⊗ M (L + n,k )) • U h -λ + n ,k (M (R + n ,k ) ⊗ M (L + n ,k )) • U h (71b)
Corollary 4.2.2. On an element interface with discontinuous but isotropic material properties on its two sides, we have [ B] = [ B ] = 0, so the characteristic coefficients {α k , α k } k=qL,qT1,qT2 of the Riemann problem (62) can be explicitly solved:

α k = M ( * L - n,k ) • (U h -U h ) , α k = M ( * L - n ,k ) • (U h -U h ) (72)
All the three fluxes defined in Theorem 4.2 are identical, and the 1st-order flux

F 1st-order n (U h , U h ) defined by (67)
is the exact solution of the Riemann problem (62).

We remark that the case of Corollary 4.2.2 was treated by Wilcox et al.

in [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF]. Our equation of numerical flux (67) is more compact and highlights the important role played by the gap between the acoustic impedances on both sides of element interfaces.

Concerning the consistency of the numerical flux (54) that is the exact solution of the Riemann problem, the following result can be proved.

Theorem 4.3. The numerical flux (69) of the first-order velocity-strain system (11) is consistent.

Proof. The second member of the system of linear equations ( 65) is equal to zero, so α n,k = α n,k = 0 for k = qT, qT 1 , qT 2 . Then, according to (64), we get

F n (U , U ) = F n (U ).

Numerical investigations of the numerical fluxes

The objectif of this section is to numerically validate the numerical fluxes proposed in the present work, especially the one solving exactly the Riemann problem and giving rise to the appropriate mechanical interface conditions ((54) in Section 3.2), which is named "MG sDG" hereafter.

In the present work, the time domain solving of the variational dG formulation (9b) is done by using the standard four-stage fourth-order Runge-Kutta iterative method, which is only conditionally stable. According to [START_REF] Cockburn | Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems[END_REF], the following stability condition written in the general context of anisotropic and piecewise homogeneous media should be satisfied:

∆t sDG ≤ CF L sDG 2 N p + 1 min E { h E max n {c E n,L } } (73
g(x, t) = a g (x) 1 -2 2π(t -T r /2) T r 2 e - 2π(t -T r /2) T r 2 (74)
with T r the period of the ricker signal. The advantage of using a ricker signal is to have a perfectly controlled frequency range with a central frequency equal

to f max = 2T -1 r
and a cutoff frequency that can be reasonably taken equal to

f c = 2.5f max .
In the present work, the highest frequency of interest is chosen equal to f c , and the associated shortest wavelength of quasi transverse waves is used to define the element size h E .

Numerical results obtained with the use of the flux "MG sDG" are compared to the reference numerical solutions calculated by using a time discontinuous space-time Galerkin solver (named "tDG" hereafter). The basic idea of this time discontinuous space-time Galerkin method is to subdivide the studied space-time domain Ω×]0, T [ into a series of space-time slabs Ω×]t n , t n+1 [ and to write in each space-time slab a variational formulation by considering both displacement and velocity unknowns. Within each space-time slab, continuous finite elements are used, but between two successive slabs, both displacement and velocity fields are discontinuous [START_REF] Hughes | Space-time finite element methods for elastodynamics: Formulations and error estimates[END_REF][START_REF] Johnson | Discontinuous Galerkin finite element methods for second order hyperbolic problems[END_REF][START_REF] Li | Implementation and adaptivity of space-time finite element method for structural dynamics[END_REF][START_REF] Tie | Adaptive time discontinuous Galerkin method for numerical modeling of wave propagation in shell and 3D structures[END_REF]. It can be shown that the time dG method is unconditionally stable, which constitutes a real advantage. Otherwise, due to energies dissipated in time jumps between two subsequent space-time slabs, there is also a numerical damping inherent to this method that increases with the frequency and allows filtering numerical spurious noises. The time dG solver used in this work discretizes the space with a finite element mesh combined with one linear element in time for each space-time slab. It can be formulated in the form of an implicit time-stepping schemes and has been validated by our previous studies [START_REF] Leclère | Parallel and adaptive FE modeling of elastic wave propagation in structures (Modélisation parallèle de la propagation d'ondes dans les structures par éléments finis adaptatifs)[END_REF][START_REF] Grédé | Numerical modeling of pyrotechnic shock wave propagation in Ariane5's structures (Modélisation des chocs d'origine pyrotechnique dans les structures d'Ariane5 : Développement de modèles de propagation et d'outils de modélisation)[END_REF][START_REF] Tie | Adaptive time discontinuous Galerkin method for numerical modeling of wave propagation in shell and 3D structures[END_REF] and also by other authors [START_REF] Li | Implementation and adaptivity of space-time finite element method for structural dynamics[END_REF][START_REF] Ekevid | Adaptive FEA of wave propagation induced by high-speed trains[END_REF].

As for the definition of the time step for the time dG solver, it is necessary to take a sufficiently small time step to prevent higher frequency modes of interest from its numerical damping. The following formula is systematically applied in our work:

∆t tDG = min E { h E max n {c E n,L } } (75)
Therefore, for all the numerical examples presented hereafter, we have ∆t tDG = 5∆t sDG , when a same element size h E is used.

Heterogeneous bimaterial case with anisotropic/isotropic materials

The first example was the one already calculated by Komatitsch et al.

[24] and de la Puente et al. [START_REF] De La Puente | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV. anisotropy[END_REF]. A square domain Ω =] -0.325, 0.325[×] -0.325, 0.325[m 2 composed of two materials is considered. The two materials are separated by the axis x 1 = 0: the one on the left is an anisotropic (transversely isotropic) zinc crystal, and the one on the right is isotropic. Their material properties are given in Table I. We remark that, for the transversely isotropic material, the symmetry axis is in the x 2 -direction. For the anisotropic mate- Not like in the works [START_REF] Komatitsch | Simulation of anisotropic wave propagation based upon a spectral element method[END_REF][START_REF] De La Puente | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV. anisotropy[END_REF], where the case of a point force acting in the x 2 -direction at location (-0.02, 0.)m was considered, a force acting in the x 1direction and defined on a segment between two points (-0.02125, 0.)m and (-0.01875, 0.)m is herein applied. The force as a time function is a ricker signal with a central frequency f max = 170kHz and a cutoff frequency f c = 425kHz.

Corresponding to the cutoff frequency, the shortest wavelength (of quasi transverse waves) is equal to 4.3mm. Two element sizes, h E 1 = 1mm and h E 2 = 0.5mm, are used here. For the "tDG" solver, the time steps used for these two mesh sizes are respectively ∆t tDG,1 = 200ns and ∆t tDG,2 = 100ns, while, for the "sDG" solver, they are respectively ∆t sDG,1 = 40ns and ∆t tDG,2 = 20ns. The total simulation time is T = 100µs.

Firstly, four snapshots of evolving quasi longitudinal and quai transverse wave fronts are presented in Figure 3 for a qualitative comparison with the analytically calculated fronts (Figure 2(a)) and the numerical results presented in [START_REF] Komatitsch | Simulation of anisotropic wave propagation based upon a spectral element method[END_REF][START_REF] De La Puente | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV. anisotropy[END_REF]. ), on the one hand, and to the results presented in [START_REF] Komatitsch | Simulation of anisotropic wave propagation based upon a spectral element method[END_REF][START_REF] De La Puente | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV. anisotropy[END_REF], on the other hand, shows that the space dG solver correctly solves all propagating phenomena at the physical interface and inside both materials.

To carry out a quantitative comparison between the tDG solver and the sDG solver, we consider displacements recorded at four sample points: S 1 = (-0.105, -0.08)m, S 2 = (-0.035, -0.08)m, S 3 = (-0.01, -0.08)m and S 4 = (0.105, -0.08)m, as proposed in [START_REF] De La Puente | An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes -IV. anisotropy[END_REF]. In Figure 3, the location of the middle point of the loading segment is indicated by a black circle, and the four sampling locations are indicated by white circles. Figure 4 show the comparison of displacement components calculated by the tDG solver and the sDG solver on the finer mesh, i.e., with h E 2 = 0.5mm. The comparison is good especially for the (quasi) longitudinal wavefronts, which is completely normal, because the wavelength of the longitudinal waves is approximately twice as long as that of the transverse waves. When the comparison is made between the results obtained by using the two different mesh sizes h E 1 and h E 2 (Figure 5), it is worth noticing that the behavior of the sDG solver is much better than that of the tDG solver. Indeed, by the tDG solver, the relative error, which is the difference between the two solutions obtained with respectively h E 1 and h E 2 , is larger, because it is mainly in the phase of the solution. In other words, the arrival time of the waves is much more correctly calculated by the sDG solver, even with the coarser element size h E 1 .

Polycrystalline materials

The second example considers the application of the sDG solver to singlephase and untextured polycrystalline materials in the 2D case and was proposed in our previous work [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF]. We recall that a single-phased polycrystal is an assembly of small single crystals bonded together. Each crystallite (grain) is anisotropic and they all have the same elastic Hooke tensor but with differently and randomly oriented basis of anisotropy. Hence, such a polycrystalline material is a typical example of piecewise homogeneous media including a large number of physical interfaces. Elastic waves propagating in it are scattered and attenuated due to interactions between waves and grain interfaces.

In a single-phase polycrystalline material, the degree of inhomogeneity in elastic properties is completely determined by the degree of anisotropy of the crystallite's Hooke tensor: the stronger the anisotropy, the higher the degree of inhomogeneity. In the related literature, it is usually expressed by measuring its departure from its Voigt average homogeneous equivalent medium in the following way [START_REF] Stanke | A unified theory for elastic wave-propagation in polycrystalline materials[END_REF]:

ξ 2 ijkl = 1 4 < (C ijkl (θ) -C V oigt ijkl ) 2 > θ (C ref -V oigt ijkl ) 2 (77)
where, the subscript "ijkl" has the proper value for either longitudinal or transverse waves, and < • > θ stands for the average over all orientations of the local material anisotropic axes denoted by θ. Each finite element model defined herein is a polycrystal occupying a rectangular domain, and it is composed of 2106 elliptic grains of size 480µm × 240µm.

The dimensions of the rectangular domain are given in Figure 6(a). The free boundary condition is prescribed on the boundary of the rectangular domain, except that a pressure loading is applied on a segment L e of length 2.8mm.

For the amplitude of the pressure loading, a Gaussian distribution along the segment L e is chosen in order to reduce as much as possible the generation of transverse waves at the ends of the segment L e . The period of ricker signal in time of the pressure loadings is T r = 0.4µs, resulting in f max = 5MHz and f c = 12.5MHz. Hence, the frequency domain of validation to be considered hereafter ranges from 2MHz to 12.5MHz.

Two polycrystals with different degrees of inhomogeneity are considered:

(1) Reference material denoted by "ref". It is a titanium alloy with crystallites of cubic symmetry. It is orthotropic and its Hooke tensor is denoted by

C ref .
(2) Material denoted by "strong". To define an orthotropic material with a high As for the choice of numerical parameters, square Q4 elements of size h = 25µm are used. This choice leads to have at least about 19 elements in the shortest quasi-longitudinal wavelength for the cut-off frequency f c for all the studied materials, except the material "strong-Voigt" for which this number is about 12. As the principle results considered are numerical measures of the attenuation coefficient and of the noise levels of the quasi longitudinal waves, the finite element size is defined with respect to the shortest quasi longitudinal wavelength, which is approximately twice the shortest quasi transverse wavelength. We note that the influence of the finite element size on these numerical measures has been studied in [START_REF] Bai | Finite element modeling of ultrasonic wave propagation in polycrystalline materials[END_REF][START_REF] Bai | Finite element modeling of grain size effects on the ultrasonic microstructural noise backscattering in polycrystalline materials[END_REF]. Finally, the time steps are calculated by using the formulas (73) for the sDG solver and (75) for the tDG solver.

C iiii C iijj,i =j C ijij,i =j ξ L C ref (
Numerical results obtained with the use of the flux "MG sDG" are compared to those calculated by either the tDG solver or the space dG solver with a flux defined by using the Voigt average of Hooke elastic tensor at a physical interface (named "cV sDG" hereafter). We recall that the flux "cV sDG" was proposed in [START_REF] Tie | A unified variational framework for the space discontinuous Galerkin method for elastic wave propagation in anisotropic and piecewise homogeneous media[END_REF] and it does not solve exactly the Riemann problem on a physical interface.

The quantitative comparison of numerical results obtained by the three solvers "tDG", "cV sDG" and "MG sDG" is made in terms of two numerically evaluated quantities: the attenuation coefficient α(f ) in the frequency domain and the backscattered noise levels in the time domain calculated as normalized root-mean-square (rms) noise levels N rms (f 0 ; t) for a given frequency f 0 . Both quantities are calculated for the longitudinal waves.

Firstly, the attenuation coefficient is numerically evaluated in the ten samples in the polycrystal "ref". In a polycrystal, the attenuation measures the amplitude decay of elastic waves during their propagation and it can be caused by dissipation, geometrical spreading or scattering-induced diffusion. In the present work, only the scattering-induced attenuation is considered, and it is usually quantified by a scalar α, called the attenuation coefficient, which defines an exponential decay law. Numerical calculation of α in the frequency domain is done in the following way [START_REF] Bai | Finite element modeling of ultrasonic wave propagation in polycrystalline materials[END_REF][START_REF] Bai | Comparison of ultrasonic attenuation within two-and three-dimensional polycrystalline media[END_REF]: Discrete Fourier Transform is used to decompose the time-series signal of the reflected wave fronts v r y (x j , t) at the jth probe x j into the frequency domain and gives rise to the corresponding amplitude spectrum vr y (x j , f ). With those data, the attenuation coefficient α(f ) as a function of the frequency f is measured as:

α(f ) = - 10 D ln( M j=1 |v r y (x j , f )| 2 M j=1 |v r,ref y (x j , f )| 2 ) ( 78
)
To eliminate the attenuation due to the geometrical spreading of wavefront, the ratio in (78) is calculated with respect to vr,ref Secondly, backscattered noise levels are numerically measured in the time domain for two frequencies f 0 = 5MHz and 10MHz. We recall that the backscat-tered noise levels are numerically measured in the time domain for a given frequency f 0 in terms of normalized root-mean-square (rms) noise levels N rms (f 0 ; t), which are defined in the following way [START_REF] Margetan | Backscattered microstructural noise in ultrasonic toneburst inspections[END_REF][START_REF] Bai | Finite element modeling of ultrasonic wave propagation in polycrystalline materials[END_REF][START_REF] Bai | Finite element modeling of grain size effects on the ultrasonic microstructural noise backscattering in polycrystalline materials[END_REF]:

N rms (f 0 ; t) =
where 10 × M means the average is taken over all probes and over all the ten samples. For a given frequency f 0 , Equation (79) defines in the time domain the rms positional average of the difference between the noise signal v b y (f 0 ; x j , t)

and the mean noise level b(f 0 ; t). The noise signal v b y (f 0 ; x j , t) is calculated by applying a frequency filtering to the noise signal v b y (x j , t) = v y (x j , t)-v ref y (x j , t)

recorded at the jth receiver x j . The normalization by E max (f 0 ) eliminates the dependence of the noise level on the incident power. E max (f 0 ) is taken equal to one half of the peak-to-peak amplitude of the incident signals for the frequency f 0 , as proposed in [START_REF] Margetan | Backscattered microstructural noise in ultrasonic toneburst inspections[END_REF].

Comparisons between the three solvers of the normalized rms noise levels N rms (f 0 ; t) in the time domain for the ten samples of the polycrystal "ref"

(resp. "strong") are presented in Figure 8 (resp. Figure 9). Backscattered noise levels evaluated by the two solvers "cV sDG" and "MG sDG" are very close.

However, it can be noted that discrepancies between them increase with time and are more due to a phase shift than to differences in the amplitude. Finally, it is shown in Figure 10 that the higher the degree of inhomogeneity, the greater the gap resulting from the use of two different numerical fluxes "cV sDG" and "MG sDG". As an important remark, the analysis of stability of a space dG method is essential, because of the use of numerical fluxes. A stability analysis has been performed by the author for the velocity-stress wave formulation [START_REF] Tie | Some comparisons and analyses of time or space discontinuous Galerkin methods applied to elastic wave propagation in anisotropic and heterogeneous media[END_REF]. At the continuous level of the time and whatever the space dimension, the stability has been proved in the case of an isotropic or anisotropic elastic medium with continuous properties or an isotropic elastic medium with discontinuous properties. But, the stability has been only proved in the 2D case for an anisotropic elastic medium with discontinuous properties satisfying a sufficient condition, which involves the degree of inhomogeneity of physical interfaces. In our current numerical studies, we have not yet encountered any instability problem. But, it is clear that this is an important point to be addressed by future studies.

  d×dsym and defined over the open set Ω×]0, T [, with R d×dsym indicating that σ is a d × d symmetric second-order tensor. The partial derivative operator with respect to time is denoted ∂ t . The space derivative

  sponding right and left eigenvectors of A n . The left eigenvectors are the eigenvectors of the adjoint of A n . As a matter of fact, among the m = d + d(d + 1)/2 eigenvalues of A n , there are d strictly negative eigenvalues (λ - n,k = -c n,k ) k=qL,{qTα} α=1,••• ,d-1 and d strictly positive eigenvalues (λ + n,k = c n,k ) k=qL,{qTα} α=1,••• ,d-1 , c n,qL and c n,qTα

Figure 1 :

 1 Figure 1: Sketch illustrating the Rankine-Hugoniot jump conditions in the Riemann problem in the 3D case.

  , [I d ] is the d × d identity matrix, and [B] and [B ] are d × d matrices with zero diagonal terms and the following extra-diagonal terms:

  ) denoted hereafter by [R] will be decomposed into two parts: [R] = [I 2d ] + [∆R]. Then the following upwind numerical fluxes are proposed: Theorem 3.2. According to how the inverse of [R] the matrix of the system of (26) is approached, three upwind numerical fluxes can be defined as follows: 1) 1st-order upwind numerical flux by taking [R] -1 ≈ [I 2d ] and by remarking that [R].[I 2d ] -[I 2d ] = [∆R]

  33) with [D], [D ], [H] and [H ] four d × d matrices defined by the following decomposition of the inverse matrix of [R]:

  49) In (49), [I d ] is the d × d identity matrix and [ B] and [ B ] are d × d matrices with zero diagonal terms and the following extra-diagonal terms:

Theorem 3 . 5 .

 35 ) as [ R] = [I 2d ] + [∆ R], the following upwind numerical fluxes are proposed: According to how the inverse of [ R] the matrix of the system of (49) is approached, three upwind numerical fluxes can be defined as follows: 1) 1st-order upwind numerical flux by taking [ R] -1 ≈ [I 2d ] and by remarking that

  ) with [ D], [ D ], [ H] and [ H ] four d × d matrices defined by the following decomposition of the inverse matrix of [ R]:[ R] -1 = [I 2d ] + [ D] -[ H] -[ H ] [ D ](55)Proof. It is straightforward by using the decomposition (15) of A n and the equations given by (46), and by substituting the solution of {α k } in the equations of numerical fluxes (47) and (48).

Corollary 3 . 5 . 1 .

 351 On an element interface with continuous material properties, we have [ B] = [ B ] = 0, ˜Ln,k = L - n,k and ˜Ln ,k = L - n ,k , so the characteristic coefficients {α k , α k } k=qL,qT1,qT2 of the Riemann problem (43) can be explicitly solved:

  the m = d + d(d + 1)/2 eigenvalues of the eigensystem (58), there are d strictly negative eigenvalues (λ - n,k = -c n,k ) k=qL,{qTα} α=1,••• ,d-1 and d strictly positive eigenvalues (λ + n,k = c n,k ) k=qL,{qTα} α=1,••• ,d-1 , c n,qL and c n,qTα being respectively the velocity of quasi longitudinal and quasi transverse wave modes propagating in the n direction. The right and left eigenmodes corresponding to these nonzero eigenvalues are: ∀k

  65) In (65), [I d ] is the d × d identity matrix and [ B] and [ B ] are the two d × d matrices already defined in Theorem 3.4. The operators M and M are defined in (12). { * L - n,k , * L - n ,k } are the perturbed left eigenmodes of the eigensystem (58) calculated by using the material properties of the adjacent element in the following way: *

Theorem 3 .Theorem 4 . 2 .

 342 4 and (65) of Theorem 4.1 have the same 2d × 2d matrix [ R], and on the other hand, the equation (66) has the same form as the equation (51) but applies to the left eigenmodes defined by (59). As in the preceding section, by decomposing [ R] into [I 2d ] + [∆ R], the following upwind numerical fluxes are proposed: According to how the inverse of [ R] the matrix of the system of (65) is approached, three upwind numerical fluxes can be defined as follows: 1) 1st-order upwind numerical flux by taking [ R] -1 ≈ [I 2d ] and by remarking that

  and [ H ] four d × d matrices already defined in Theorem 3.4. Proof. It is straightforward by using the decomposition (61) of A n and by substituting the solution of {α k } in the equation of numerical flux (64).

  ) with N p the order of the FE basis function, h E the size of the element E and CF L sDG the Courant-Friedrichs-Levy number. For the numerical examples presented hereafter, FE meshes with four-node quadrilateral (Q4) elements are systematically used. Hence, we have N p = 1. The constant CF L sDG depends on the order of FE basis function, the space dimension d (d = 1, 2, 3), as well as the shape of finite element. According to our numerical experiences, CF L sDG = 0.6 always works in the 2D case with N p = 1, i.e., for three-node triangular (T3) elements or Q4 elements. The element size h E should be chosen by considering the shortest wavelength of interest, which is, in turn, defined by the highest frequency of interest often determined by the frequency content of external loadings. For the numerical examples presented hereafter, ricker signals are used to define the dependency of the external loadings upon the time t, i.e.:

Figure 2 :

 2 Figure 2: (a) Wave velocities, and (b) Waves fronts, of the anisotropic material (only curves in red blue are to be considered in the 2D problem defined in the Ox 1 x 2 plan).

Figure 3 :

 3 Figure 3: (a) and (b) divu respectively at 30µs and 60µs; (c) and (d) curlu respectively at 30µs and 60µs. The location of the middle point of the loading segment is indicated by a black circle, and the four sampling locations are indicated by white circles.

Figure 4 :

 4 Figure 4: Comparison of displacement components calculated by the tDG solver and the sDG solver on the finer mesh.

Figure 5 :

 5 Figure 5: Comparison of displacement components calculated with two different mesh sizes.

Figure 6 :

 6 Figure 6: (a) FE model of a single-phase polycrystal composed of elliptic grains with pressure loading applied on the emitter segment Le with a Gaussian distribution. Dispersions dV /Vm in the phase wave velocity c qL (k) with k = ey for the polycrystals (a) "ref" and (b) "strong", Vm being the averaged wave velocity < c qL (k) > grains over all grains and dV = c qL (k) -Vm for each grain.

y

  (x j , f ) the reflected echo signals recorded in the equivalent homogeneous medium instead of the incident signals.In (78), M = 22, the constant 10 is for the unit conversion from Neper (N p) to decibel (dB), and D stands for the wave propagation distance just before the arrival of reflected echoes at the probes.

Figure 7 9 (

 79 Figure 7 presents comparisons between the numerical measurements of attenuation in the polycrystal "ref" by the three solvers 'tDG", "cV sDG" and "MG sDG". The attenuation coefficient α is plotted either for the crystallographic orientation distribution Θ 0 (Figure 7(a)) or for all the ten {Θ i } i=0,••• ,9(Figure7(b)). In the latter case, the averaged measure and the minimum and maximum bounds for all the ten measures are plotted. The numerical measurements given by both space dG solvers with two different fluxes are very close, giving a maximum relative difference of 3% for both single and averaged measurements. We note that the polycrystal "strong" is not considered for the attenuation coefficient, because coherent echo signals cannot be recorded due to the high levels of backscattered noises.

Figure 7 :

 7 Figure 7: Comparison between three solvers of numerical measurements of the attenuation coefficient α in the polycrystal "ref". (a) α measured in the sample Θ 0 ; (b) Discrepancies of numerical measure with mean Θ (α), min Θ (α) and max Θ (α) denoting respectively the averaged numerical measures, the minimum and maximum bounds of numerical measures over the ten samples Θ 0 -Θ 9 ; (c) Relative differences taking the measurements given by "MG sDG" as the reference one.

  b y (f 0 ; x j , t) -b(f 0 ; t)) 2 E max (f 0 ) 0 ; x j , t)

Figure 8 :

 8 Figure 8: Comparison between three solvers of the normalized rms noise levels Nrms(f 0 ; t) for the ten samples of the polycrystal "ref". (a) f 0 = 5MHz; (b) f 0 = 10MHz

Figure 9 :

 9 Figure 9: Comparison between three solvers of the normalized rms noise levels Nrms(f 0 ; t) for the ten samples of the polycrystal "strong". (a) f 0 = 5MHz; (b) f 0 = 10MHz

Figure 10 :

 10 Figure 10: Comparison between both solvers "cV sDG" and "MG sDG" of the averaged normalized rms noise levels Nrms(f 0 ; t): relative differences taking the measurements given by "MG sDG" as the reference one for the polycrystals (a) "ref" and (b) "strong"

Table II :

 II Elastic moduli and degrees of inhomogeneity of the studied materials "ref" and "strong", and their homogenized equivalent media "ref-Voigt" and "strong-Voigt"

	GPa)	134.0	110.0	36.0	0.028
	C ref -V oigt (GPa) 153.0	100.0	26.5	0
		C 1111 C iiii,i=2,3 C 2233 C 11ii,i=2,3 C 1i1i,i=2,3 C 2323	ξ L
	C strong (GPa)	132.5	10.3	0.4	3.6	4.0	5.0	0.563
	C strong-V oigt (GPa)	35.4	35.4	10.6	10.6	12.4	12.4	0

Acknowledgments

This work was granted by the ANR (Agence Nationale de la Recherche) -MAPIE project, ANR-13-MONU-002. Computations were performed using HPC resources from the "Mésocentre" computing center of CentraleSupélec and ENS Paris-Saclay.

Appendix A1

Proof of Theorem 3.1. By adding the three equations of [START_REF] Leclère | Parallel and adaptive FE modeling of elastic wave propagation in structures (Modélisation parallèle de la propagation d'ondes dans les structures par éléments finis adaptatifs)[END_REF], i.e. Eq.(21a) + Eq.(21b) + Eq.(21c), we get:

According to the definition of the Jacobian operator [START_REF] Wilcox | A high-order discontinuous Galekin method for wave propagation through coupled elasticacoustic media[END_REF], the eigenmodes [START_REF] Cockburn | Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems[END_REF] and the definition [START_REF] Johnson | Discontinuous Galerkin finite element methods for second order hyperbolic problems[END_REF] and the property [START_REF] Hughes | Space-time finite element methods for elastodynamics: Formulations and error estimates[END_REF] of the Christoffel tensor Γ n , we obtain the following equations by considering separately ℘ vect (Eq.( 81)) and n • ℘ tens (Eq.( 81)):

We recall that the sum is taken over l = qL, qT 1 , qT 2 and for the sake of clarity the symbol of summation l is explicitly added in (82). ( 82) is a linear system of six equations for six unknowns {α k , α k }. Then, to obtain the expressions given in Theorem 3.1, the following manipulations of (82) are made:

• For first three equations corresponding to

• For the last three equations corresponding to

The matrix form of ( 83) and (84) gives rise to [START_REF] Bai | Finite element modeling of ultrasonic wave propagation in polycrystalline materials[END_REF].

Appendix A2

Lemma A2.1. In the case of continuous material properties, the following equations hold on the interface of elements E and E :

Proof. As the material properties are continuous, we have Γ n = Γ n , so λ n ,k = λ n,k . Otherwise, by recalling ( 16), it can be verified that:

k can be verified in the same way.

Appendix A3

Proof of Theorem 3.4. By adding the three equations of (43), i.e. Eq.(43a) + Eq.(43b) + Eq.(43c), we get:

According to the definition of Ãn (45), the eigenmodes ( 16) and the equation (46), we obtain the following equations by considering separately ℘ vect (Eq.( 87)) and ℘ tens (Eq.( 87)):

To obtain (88b), we have used the fact that w ⊗ s n = 0 implies w = 0.

Then, to obtain the expressions given in Theorem 3.4, the following manipulations of (88) are made:

• For the first three equations corresponding to

Eq.(88a) + Eq.( 88b)) leads to:

• For the last three equations corresponding to

Eq.(88a) + Eq.(88b)) leads to:

The matrix form of ( 89) and (90) gives rise to (49).

Appendix A4

Proof of Theorem 4.1. By adding the three equations of (62), i.e. Eq.(62a) + Eq.(62b) + Eq.(62c), we get:

According to the definition of A n (13) and the eigenmodes [START_REF] Cockburn | Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems[END_REF], we obtain the following equations by considering separately ℘ vect (Eq.( 91)) and ℘ tens (Eq.( 91)):

As for (88b), the fact that w ⊗ s n = 0 implies w = 0 is used to obtain (92b).

It is interesting to remark that (92) is identical to (88) if we change C : ε h with σ h .

Then, to obtain the expressions given in Theorem 4.1, the following manipulations of (92) are made:

• For the first three equations corresponding to

Eq.(92a) + Eq.(92b)) leads to: Eq.(92a) + Eq.(92b)) leads to:

The matrix form of ( 93) and (94) gives rise to (65).