Yu Liu

Benhui Fan

Ahmed Sridi

Angkeara Svay

Khalil Abene

Luciano De

Carvalho Paludo

Lucio De

Abreu Correa

Fillipo Gatti

Maroua Hammani

Alfonso Panunzio

Hélène Moustacas

Matthieu Bonneric

Delong He

Milanetto Thiago

Romain Schlittler

Keywords:

First and foremost, I would like to express my sincere gratitude to my supervisors Bing Tie and Régis Cottereau, and also Denis Aubry

Context

Wave propagation in homogeneous media and its numerical modelling have been widely developed for many years [START_REF] Achenbach | Wave propagation in elastic solids[END_REF][START_REF] Bedford | Elastic wave propagation[END_REF][START_REF] Karal | Elastic wave propagation in homogeneous and inhomogeneous media[END_REF][START_REF] Mccarthy | Elastic Wave Propagation[END_REF]. However, when waves propagate in heterogeneous media, where the physical properties of media vary continuously or discontinuously, waves interact with heterogeneities and a lot of complicated phenomena arise: reflection, refraction, diffraction, scattering, interference, attenuation, dispersion, etc. [START_REF] Achenbach | Wave propagation in elastic solids[END_REF][START_REF] Ishimaru | Wave propagation and scattering in random media[END_REF][START_REF] Manolis | Seismic Wave Propagation in Non-Homogeneous Elastic Media by Boundary Elements[END_REF][START_REF] White | Seismic waves: Radiation, transmission, and attenuation[END_REF]. Due to this complex nature, direct numerical simulations on the fine grids become prohibitively expensive and the quantities that are classically considered as relevant for quality control of numerical solutions are always highly oscillating in time-space. It is therefore important and necessary to develop efficient, accurate and even new methods for the quantification of numerical errors. This work deals with the numerical analysis of wave propagation problems in heterogeneous media and the a posteriori error estimation for its finite element solutions, which allows assessing and controlling their quality. This project has various applications in different domains, such as the numerical modeling of polycrystalline permeability to ultrasonic waves [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][START_REF] Hirsekorn | The scattering of ultrasonic waves by polycrystals[END_REF][START_REF] Van Pamel | Finite element modelling of elastic wave scattering within a polycrystalline material in two and three dimensions[END_REF]. Crystalline solids are composed of large quantities of grains (or crystallites). Each grain can belong to a specific class of crystal symmetry exhibiting anisotropic elasticity and possess a unique crystallographic orientation from the neighbouring grains. An elastic wave travelling through a polycrystalline material gets scattered at grain interfaces due to crystallographic misorientations of adjacent grains. This scattering leads to a loss of energy by the coherent wave and results in apparent wave attenuation and noise signals, which depend on ultrasonic frequency, grain size and shape, and the elastic properties of the crystallites. Understanding elastic wave interactions with these media is very important for non-destructive testing [START_REF] Blitz | Ultrasonic methods of non-destructive testing[END_REF]. Figure 1.1 illustrates ultrasonic wave propagation in a polycrystalline material.

This project can also be applied in geophysics [START_REF] Dhua | Wave propagation in heterogeneous layers of the earth[END_REF][START_REF] Sato | Seismic wave propagation and scattering in the heterogeneous earth[END_REF][START_REF] Virieux | Modelling seismic wave propagation for geophysical imaging[END_REF]. Seismic or elastic wave propagation through the earth plays an important role in understanding damages during earthquakes. It is known that the materials in the earth are highly heterogeneous and their elastic properties vary with the depth from one region to another. Although this variation may be gradual, there exist also discontinuities that separate media with different elastic coefficients. The existence of arbitrary-shaped interfaces, inclusions and fractures requires an efficient treatment [START_REF] Aggelis | Numerical simulation of surface wave propagation in material with inhomogeneity: Inclusion size effect[END_REF][START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF][START_REF] Lombard | Modeling 1-d elastic p-waves in a fractured rock with hyperbolic jump conditions[END_REF]. Considering that the earth has strong heterogeneities and complex geometries, its numerical modelling is quite a challenging problem. Figure 1.2b gives P-wave velocity of the Marmousi2 model. The geometry of this model is based on a profile through the North Quenguela trough in the Cuanza basin [START_REF] Brougois | Marmousi, model and data[END_REF][START_REF] Martin | Marmousi2: An elastic upgrade for marmousi[END_REF]. Figure 1.2c shows a simulation of P-wave propagation in soil. Furthermore, in geophysics, the length scale is of the order of the kilometre and the physical properties vary continuously or discontinuousely, while in the polycrystalline application, the length scale is of the order of the micrometre and the anisotropic material properties are discontinuous from one grain to another. Even if the involved scales and the frequency of waves are very different in the two cases and the polycrystalline materials reveal much more discontinuities at the interfaces, they present similar issues when efficient numerical simulations are sought for.

Modelling of elastic wave propagation in heterogeneous media

In this section, we start with the introduction of the elastodynamic equation in homogeneous media. Then different types of waves and various propagation phenomena in heterogeneous media are presented. Finally, the high frequency (HF) wave propagation in the weak coupling regime studied in this work is introduced in detail.

Wave propagation in homogeneous media

We introduce firstly the elastic wave propagation in homogeneous media in an open domain Ω ⊆ R 3 and R 3 stands for the general three-dimensional Euclidean space. The density of media is denoted by ρ. The displacement field is denoted by u(t, x) : (0, T ) × Ω → R 3 , where x ∈ Ω is the spatial position vector. The second-order Cauchy stress tensor is denoted by σ σ σ (t, x) and f stands for the density of body force.

Then the second-order elastodynamic equation reads:

ρ∂ 2 t u(t, x) -∇ x • σ σ σ (u(t, x)) = f ∀(t, x) ∈ (0, T ) × Ω (1.1)
where ρ is independent of x in the homogeneous case and ∇ x is the gradient with respect to x. The stress field σ σ σ is defined by the generalized Hooke's law:

σ σ σ (t, x) = C : ε ε ε(t, x) (1.2)
where C is the fourth-order elasticity tensor, independent of x in the homogeneous case. For homogeneous isotropic elastic materials [START_REF] Achenbach | Wave propagation in elastic solids[END_REF], eq. (1.2) becomes:

σ σ σ = λ tr(ε ε ε)I 3 + 2µε ε ε (1.3)
where tr(•) is the trace operator, I 3 is the 3 × 3 identity matrix, and λ , µ are two Lamé parameters.The strain tensor ε ε ε is defined by:

ε ε ε(t, x) = ∇ x ⊗ s u(t, x) = 1 2 ∇ x ⊗ u(t, x) + (∇ x ⊗ u(t, x)) T (1.4)
where ⊗ s is the symmetrized tensor product of two vectors, ⊗ is the tensor product and the superscript T stands for the transpose.

Equation (1.1) is solved subject to the initial conditions and the boundary conditions. At initial instant t = 0, the displacement and the velocity are defined as:

u(0, x) = u 0 (x), ∂ t u(0, x) = v 0 (x) ∀x ∈ Ω (1.5)
1.2 Modelling of elastic wave propagation in heterogeneous media 7

The boundary conditions mean that the displacement and/or the stress vectors are specified on ∂ Ω, such as Dirichlet or Neumann boundary conditions [START_REF] Weisstein | Boundary conditions. From MathWorld-A Wolfram Web Resource[END_REF]:

u = 0 on (0, T ) × Γ u , σ σ σ (u) • n = h on (0, T ) × Γ σ σ σ (1.6)
where h denote the density of surface force.

In the particular case of harmonic plane waves, the displacement fields have the following form:

u(t, x) = Ue i(x•k-tω) (1.7)
where U is the wave mode defining the polarization direction and amplitude, i is the imaginary unit, k is the wave vector which specifies the direction of propagation and ω is the angular frequency of wave. Now introducing eq. (1.7) to eq. (1.1) and taking f = 0, the so-called Christoffel equation is obtained:

ω 2 I 3 -Γ Γ Γ(k) • U = 0 (1.8)
where the Christoffel tensor Γ Γ Γ is defined as (with Einstein summation convention): ∀W,

Γ Γ Γ(k) • W = ρ -1 (C : (W ⊗ s k)) • k i.e. Γ ik = ρ -1 C i jkl k j k l (1.9)
From eq. (1.8) it is known that the polarization vector U and ω 2 are respectively the eigenvector and eigenvalue of Γ Γ Γ(k). The wave propagation modes can be obtained by using eigendecomposition. For homogeneous isotropic elastic materials, Γ Γ Γ(k) has three eigenvalues:

ω 2 1 = λ + 2µ ρ |k| 2 , ω 2 2 = ω 2 3 = µ ρ |k| 2 (1.10)
The corresponding eigenvectors or polarization directions of waves are respectively denoted by three unit vectors k, k⊥ 1 and k⊥ 2 . k = k |k| is the eigenvector of ω 2 1 and k⊥ 1 , k⊥ 2 are the eigenvectors of ω 2 2 such as ( k, k⊥ 1 , k⊥ 2 ) form an orthonormal triplet. Equation (1.10) are called the dispersion relation for each wave mode. The first eigenvalue corresponds to the P-wave mode and the other one with multiplicity of two to the S-wave mode:

ω p = c p |k| = λ + 2µ ρ |k| , ω s = c s |k| = µ ρ |k| (1.11)
where c p and c s are respectively velocities of P-wave and S-wave. P-waves propagate at a higher velocity than do the S-waves and the two S-waves in an isotropic medium propagate at the same velocity. surface [START_REF] Achenbach | Wave propagation in elastic solids[END_REF]151]. There are two types of surface waves: Love waves and Rayleigh waves [START_REF] Aki | Quantative seismology: Theory and methods[END_REF]. The Love wave is polarized in the horizontal direction perpendicular to the propagation path, while Rayleigh waves have an ellipsoidal polarization in the vertical plane through the path of propagation (fig. 1.3, right).

Phenomena and scattering regimes for wave propagation in heterogeneous media

Now we focus on heterogeneous media, where physical properties of media C depend on x (ρ can also depend on x but here we assume it constant).

As we mentioned in the section 1.1, various complex phenomena appear when different types of waves propagate in heterogeneous media. More specifically, waves are reflected when they encounter boundaries of material properties through which they are traveling. When waves meet the interface between two media of different material properties, they change directions and this phenomenon is called refraction. Diffraction refers to the ability of waves to bend around corners or slits and to spread whenever they encounter obstacles. Attenuation is the term used to account for loss of wave amplitude due to all mechanisms, including absorption, scattering, and mode conversion. Scattering refers to wave radiation from heterogeneity acting as secondary sources of radiation due to excitation by the incident wave. An example of illustration of some phenomena mentioned above is given in fig. 1.4. Specifically, scattering can be classified into different regimes according to the characteristic length scales [START_REF] Aki | Scattering and attenuation of seismic waves, part i[END_REF][START_REF] Wu | Introduction: Seismic wave scattering in three-dimensionally heterogeneous earth[END_REF]. In fact, it is known that wave propagation in heterogeneous media is a multi-scale problem both mathematically and numerically, and different phenomena can happen for different scales involved. There are at least three fundamental length scales for wave propagation problems: the propagation distance L, the characteristic length of heterogeneity l c (the scale on Fig. 1.4 Acoustic wave propagation phenomena: specular reflection, diffuse reflection, refraction and diffraction [START_REF] Kapralos | Sonel mapping: a probabilistic acoustical modeling method[END_REF] which the heterogeneous medium varies) and the dominant wavelength λ . Also, the amplitude of the fluctuations of heterogeneities σ 2 has to be considered. The following scattering regimes and phenomena can be described with different values of dimensionless parameters ξ := kl c = 2πl c /λ and ε := λ /L [START_REF] Yoon | Deep seismic imaging in the presence of a heterogeneous overburden[END_REF]:

1. quasi-homogeneous or effective medium regime: ξ < 0.01, σ 2 ≈ 1 and ε ≈ 1. The medium can be regarded as quasi-homogeneous, because the heterogeneous scale length l c is very small compared to the typical wavelength λ . The homogenization methods with a deterministic effective wave equation can be applied [START_REF] Capdeville | 1-d non-periodic homogenization for the seismic wave equation[END_REF][START_REF] Fish | Space-time multiscale model for wave propagation in heterogeneous media[END_REF][START_REF] Lombard | Numerical modeling of the acoustic wave propagation across a homogenized rigid microstructure in the time domain[END_REF]. Introduction 4. forescattering regime: ξ ≫ 1. The wavelength is small compared to the length of heterogeneity and the scattered wave energy is mainly concentrated near the forward direction. In this case the parabolic approximation for the full wave equation can be used [START_REF] Corones | A new parabolic approximation to the helmholtz equation[END_REF]. It describes a one-way (forward direction) wave propagation, since the back-scattered waves are very weak. The scattering problem becomes a focusing diffraction and interference problem [START_REF] Li | Influence of forward scattering on ultrasonic attenuation measurement[END_REF][START_REF] Li | Geometrical-optics approximation of forward scattering by gradient-index spheres[END_REF].

High frequency setting

The high frequency wave propagation, when we assume the wavelength is small compared to the propagation distance (i.e. ε ≪ 1), is specifically introduced here. It should be studied in a different and specific way compared to the low or mid frequency wave propagations. The most relevant observables for modelling wave propagation in HF limit are no longer the wave field u ε or its derivatives, but rather its quadratic quantities like energy or energy density. It can be understood by a simple example discussed in [START_REF] Baydoun | Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media[END_REF][START_REF] Savin | Kinetic modeling for transport of elastic waves in anisotropic heterogeneous media[END_REF].

Considering the following wave field x → u ε (x), which oscillates with an amplitude a(x) and its mean function u(x):

u ε (x) = u(x) + a(x)sin x ε , 0 < ε ≪ 1 (1.12)
u ε (x) has no strong limit when ε → 0, although a(x) and u(x) vary slowly. However, for any continuous function ϕ(x) having a compact support on R, a weak limit can be obtained:

lim ε→0 R ϕ(x)(u ε (x)) 2 dx = R ϕ(x) u 2 (x) + 1 2 a 2 (x) dx (1.13)
as a consequence of Riemann-Lebesgue's lemma. The observation function ϕ(x) allows computing a smoothed "energy" of u ε (x), locally at a selected point of interest x, given by u 2 (x) + 1 2 a 2 (x), and quantifying the influence of oscillations with amplitude a at that point (fig. 1.5). Above all, this weak limit has no longer the oscillating properties of u ε (x).

Fig. 1. [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous galerkin finite element methods for the second-order wave equation[END_REF] The function u ε (x) (thin solid line), its mean u(x) (thick solid line) and its square root weak limit (u 2 (x) + 1 2 a 2 (x))

1 2 (thick dashed line) [START_REF] Baydoun | Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media[END_REF] 1.2 Modelling of elastic wave propagation in heterogeneous media 11

Similarly, for wave propagation problems, the highly oscillatory features of quantities such as the velocity and pressure fields (for acoustic waves) or the displacement and stress fields (for elastic waves) are more difficult to understand in a high frequency regime. The energy or the energy density arise as the most relevant observables for the characterization of HF wave propagation. These quadratic quantities are no more highly oscillating. More details can be found in [START_REF] Baydoun | Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media[END_REF][START_REF] Savin | Kinetic modeling for transport of elastic waves in anisotropic heterogeneous media[END_REF][START_REF] Zeng | Modeling of high-frequency seismic-wave scattering and propagation using radiative transfer theorymodeling of high-frequency seismic-wave scattering and propagation using radiative transfer theory[END_REF]. In this context, we shall thus consider the evolution of the energy density (i.e. Wigner measure, presented in detail in the chapters 3 and 4) associated to high frequency acoustic or elastic waves.

Kinetic model for high frequency wave propagation in the weak coupling limit

When high frequency elastic waves propagate in random media in the weak coupling limit, i.e. λ ≈ l c ≪ L, σ 2 ≪ 1 (1.14) and the background physical properties of weakly heterogeneous media can be modelled by a superposition of a homogeneous or slowly varying part, called also slowly background, and a fast fluctuating part, a kinetic model which describes the evolution of wave energy will arise.

Its derivation can be based on the use of the Bethe-Salpeter equation for the propagation of correlations, which satisfies the principle of conservation of total wave energy and the second law of thermodynamics [START_REF] Howe | On the kinetic theory of wave propagation in random media[END_REF][START_REF] Howe | A kinetic equation for wave propagation in random media[END_REF]. Or more precisely mathematically, a multi-scale asymptotic analysis by Wigner transforms and their interpretation in terms of semiclassical operators can also be used [START_REF] Bal | Kinetics of scalar wave fields in random media[END_REF][START_REF] Savin | Kinetic modeling for transport of elastic waves in anisotropic heterogeneous media[END_REF]. For example, Ryzhik et al. [START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF] derived and analyzed radiative transfer equations for the energy density of acoustic, electromagnetic, and elastic waves in random media with the spatial Wigner transform; Baydoun et al. [START_REF] Baydoun | Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media[END_REF] assessed the influence of material anisotropy on the radiative transfer regime of elastic waves in randomly heterogeneous media. All these methods allow deriving a transport equation (for homogeneous media) or a radiative transfer equation (for heterogeneous media) from the wave equation in terms of energy densities in phase space. We present the main derivation steps of this equation obtained by Wigner transform for homogeneous and heterogenous media respectively in the sections 3.2.1 and 4.1.2.

In general, analytical solutions do not exist for the radiative transfer equation. However, when considering scalar wave propagation in an isotropic homogeneous background medium, it can be solved analytically [START_REF] Howe | A kinetic equation for wave propagation in random media[END_REF][START_REF] Paasschens | Solution of the time-dependent boltzmann equation[END_REF][START_REF] Van Wijk | 1d energy transport in a strongly scattering laboratory model[END_REF].

Bal and Pinaud [START_REF] Bal | Accuracy of transport models for waves in random media[END_REF] compared particularly the energy densities of wave fields by solving the acoustic wave equation with those obtained by solving the radiative transfer equation. The wave equation was simulated by a finite difference forward scheme in time and the radiative transfer equation was solved numerically by Monte Carlo method [START_REF] Lapeyre | Introduction to Monte-Carlo methods for transport and diffusion equations[END_REF][START_REF] Margerin | Monte carlo simulation of multiple scattering of elastic waves[END_REF]. Finally a good agreement was found between the two energy quantities integrated in a subdivision of the studied 2D random media. It illustrates the equivalence of these two equations in the weak coupling regime and allows validating the weak limit proposed by the theory of radiative transfer in random media.

As we presented above, the energy density is a more relevant observable in the high frequency limit than wave fields and the radiative transfer equation is equivalent to the wave equation in random media in terms of energy. Thus in this work, we aim at studying this new equation and quantifying errors of numerical solutions of the wave equation in terms of energy based on it.

Error estimation for numerical solutions of wave equation

For wave propagation in heterogeneous media with complex phenomena presented in the section 1.2.2, described by the partial differential equation (PDE) with varying coefficients and initial or boundary conditions, it is typically impossible to find exact analytical solutions. Various numerical methods are used to find an approximate numerical solution. It is important to verify their validity and quality by some tools of error estimation.

Numerical methods for wave equation

In dynamics, the usual discretization of problem combines space discretization and a time-stepping scheme. The most common numerical methods for space discretization in elastodynamics are the finite difference method (FDM) [START_REF] Boore | Finite difference methods for seismic wave propagation in heterogeneous materials[END_REF][START_REF] Smith | Numerical solution of partial differential equations: finite difference methods[END_REF][START_REF] Virieux | P-sv wave propagation in heterogeneous media: Velocity-stress finitedifference method[END_REF] and the finite element method (FEM) [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF][START_REF] Seron | Finite-element method for elastic wave propagation[END_REF][START_REF] Smith | The application of finite element analysis to body wave propagation problems[END_REF][START_REF] Zienkiewicz | The finite element method[END_REF]. Others include the spectral method [START_REF] Gottlieb | Numerical analysis of spectral methods: theory and applications[END_REF][START_REF] Taylor | Spectral methods for the wave equation in second-order form[END_REF], the boundary-element method (BEM) [START_REF] Banerjee | Boundary element methods in engineering science[END_REF][START_REF] Mansur | Numerical implementation of the boundary element method for two dimensional transient scalar wave propagation problems[END_REF], the finite-volume method (FVM) [START_REF] Dormy | Numerical simulation of elastic wave propagation using a finite volume method[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF].

More specifically, the FDM, based directly on the strong formulation of PDE, discretizes the exact solution through approximations of the partial differential operators. It is mainly adapted to simple geometries. The FEM is based on the variational or the weak form of PDE. It is better adapted to complicated geometries, but sometimes requires more memory space for numerical simulations than FDM. The spectral element method (SEM) [START_REF] Gatti | Investigation of the earthquake ground motion coherence in heterogeneous non-linear soil deposits[END_REF][START_REF] Komatitsch | Introduction to the spectral element method for threedimensional seismic wave propagation[END_REF] is a formulation of FEMs with a subdivision of the computational domain into hexahedral elements. The most important property of the SEM is that the mass matrix is diagonal by construction, which saves time and memory. The BEM transforms PDE describing a boundary value problem to an equivalent representation by integral equations with known and unknown boundary states. Hence, it only requires discretization of the boundary surface rather than the volume, i.e., the dimension of problems is reduced by one. However, the BEM matrix is a unsymmetric and full matrix with non-zero coefficients and the FEM matrix is much larger but very sparsely populated. It means that the latter can be stored and solved more efficiently. Besides, the transformation of PDE to boundary integral equations requires the use of Green function, which is difficult and restricts the applications of the BEM, especially for problems in heterogeneous media. Based on the strong formulation of PDE, the FVM divides the domain into elements (called "control volume"). By the divergence formula, an integral formulation of the fluxes over the boundary of the control volume is then obtained. The fluxes on the boundary are discretized with respect to the discrete unknowns. The FVM is usually used in solving fluid flow problems. As the FEM, the advantage of the FVM is that it is easily formulated to allow for unstructured meshes.

The FEM is chosen in this work. It has been widely used for solving wave propagation problems [START_REF] Dhia | Time-harmonic acoustic propagation in the presence of a shear flow[END_REF][START_REF] Van Pamel | Finite-element modelling of elastic wave propagation and scattering within heterogeneous media[END_REF]. Here we present in brief the discretization for the wave equation with the FEM and the related notations and definitions.

Error estimation for numerical solutions of wave equation 13

Firstly, the weak form of the wave equation (1.1) with initial conditions (1.5) and boundary conditions (1.6) is built by introducing a test function w:

∀t, find u(t, •) ∈ V(Ω) with u(0, x) = u 0 (x), ∂ t u(0, x) = v 0 (x) such that ∀w ∈ V(Ω), (∂ 2 t u, w) M(Ω) + (u, w) K(Ω) = (f, w) (Ω) + ⟨h, w⟩ (Γ σ σ σ ) (1.15)
where V(Ω) is the space of kinematically admissible displacements. The products in eq. (1.15) are defined as:

(u, w) (Ω) = Ω u • wdx, (u, w) K(Ω) = Ω σ σ σ (u) : ε ε ε(w)dx, (u, w) M(Ω) = (ρu, w) (Ω) , ⟨u, w⟩ (Γ σ σ σ ) = Γ σ σ σ u • wdx (1.16)
Let V h (Ω) be a subspace of V(Ω), the FEM consists in reconstructing solutions expanded on a basis {w h } of V h (Ω), which is defined using a spatial mesh T h , a partition of the domain Ω into elements (S). The letter h stands for the maximal diametre of the elements. A semi-discretized weak form of the variational problem (1.15) reads:

∀t, find u h (t, •) ∈ V h (Ω) with u h (0, x) = u 0,h (x), ∂ t u h (0, x) = v 0,h (x) such that ∀w h ∈ V h (Ω), (∂ 2 t u h , w h ) M(Ω) + (u h , w h ) K(Ω) = (f, w h ) (Ω) + ⟨h, w h ⟩ (Γ σ σ σ ) (1.17)
The approximate displacement field u h is defined by:

u h (t, •) = ∑ A,k U Ak (t)w A (x)e k (1.18)
where A is the node of mesh, k stands for the direction in space, e k is the unit vector in direction k, w A is the shape function at node A, U Ak is the node displacement at node A and in direction k. Equation (1.17) can be expressed in the following matrix form:

[M] Ü(t) + [K] {U(t)} = {F(t)} (1.19)
where {U(t)} is the node displacement vector for instant t, [M] is the mass matrix and [K] is the stiffness matrix defined respectively by:

M (Ak)(Bl) = (w A e k , w B e l ) M(Ω) , K (Ak)(Bl) = (w A e k , w B e l ) K(Ω) (1.20)
and {F(t)} is the vector of forces:

F (Bl) = (f, w B e l ) (Ω) + ⟨h, w B e l ⟩ (Γ σ σ σ ) (1.21)
For solving the system (1.19), various time-stepping schemes can be chosen for the second-order time derivative, such as the forward finite difference scheme [START_REF] Georgoulis | A posteriori L ∞ (L 2 )-error bounds for finite element approximations to the wave equation[END_REF], the Newmark scheme [START_REF] Marfurt | Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations[END_REF], etc.

Introduction

Take the forward finite difference scheme as an example. The interval (0, T ) is divided into equal subintervals I n = (t n ,t n+1 ) (n ∈ [0, N -1] with t 0 = 0 and t N = T ) of length ∆t. We denote u n h∆t = u h (t n ) (u 0 h∆t = u 0 ), the time-stepping scheme reads:

∂ 2 t u n h∆t =    ∂ t u n h∆t -∂ t u n-1 h∆t ∆t , n ∈ [1, N] a 0 , n = 0 with ∂ t u n h∆t =    u n h∆t -u n-1 h∆t ∆t , n ∈ [1, N] v 0 , n = 0 (1.22)
Then time discrete solutions u n h∆t are found for all n with eq. (1.22) and eq. (1.19). The semidiscretization in space and time-stepping scheme presented above can also be based on a first-order equation by introducing u = ∂ t v.

Space-time finite element methods, with or without the continuity in time, can also be used in elastodynamics. For instance, the time discontinuous space-time Galerkin (DG) method has been developped for solving the wave equation [START_REF] Hughes | Space-time finite element methods for elastodynamics: formulations and error estimates[END_REF][START_REF] Li | Structural dynamic analysis by a time-discontinuous galerkin finite element method[END_REF][START_REF] Tie | Adaptive time discontinuous galerkin method for numerical modelling of wave propagation in shell and 3d structures[END_REF]. It is based on the following strong first-order form of eq. (1.1):

   ρ∂ t v(t, x) -∇ x • σ σ σ (u(t, x)) = f ∇ x • (σ σ σ (∂ t u(t, x) -v(t, x))) = 0 (1.23)
The time period is subdivided as in what we presented for the time-stepping scheme above. Let a finite element space V n (M n ) is given for each space-time slab M n = I n × Ω. The weak form of DG method is then formulated as follows:

find (u, v) ∈ ∏ N-1 n=0 V n (M n ) × V n (M n ) such that ∀w = (w u , w v ) ∈ ∏ N-1 n=0 V n (M n ) × V n (M n ) A ((u, v), w) = L(w) (1.24) with A ((u, v), w) = N-1 ∑ n=0 I n (∂ t v, w v ) M(Ω) dt + N-1 ∑ n=0 I n (u, w v ) K(Ω) dt equilibrium + N-1 ∑ n=1 [v] n , w n v,+ M(Ω) dt + v 0 + , w 0 v,+ M(Ω) continuities of v in time + N-1 ∑ n=0 I n (∂ t u, w u ) K(Ω) dt - N-1 ∑ n=0 I n (v, w u ) K(Ω) dt equivalence of ∂ t u and v + N-1 ∑ n=1 I n ([u] n , w n u,+ ) K(Ω) dt + (u 0 + , w n u,+ ) K(Ω) continuities of u in time , L(w) = v 0 , w 0 v,+ M(Ω) + (u 0 , w 0 u,+ ) K(Ω) initial conditions + N-1 ∑ n=0 I n (f, w v ) (Ω) dt + N-1 ∑ n=0 I n ⟨h, w v ⟩ (Γ σ σ σ ) dt boundary conditions (1.25) where u n ± (x) = lim τ→0 + u(t n ± τ, x), [u] n = u n + -u n -.
According to eq. (1.24), the displacement and velocity fields (u, v) are continuous in space and discontinuous in time between two successive space-time slabs and the continuity in time is enforced only weakly. To solve eq. (1.24) numerically, the time-space domain M n is discretized by 1.3 Error estimation for numerical solutions of wave equation a space-time finite element mesh, and u n h∆t and v n h∆t are then obtained. In our studies, the choice of V n (M n ) = P 1 (I n ) × V n (Ω) is adopted, where P 1 (I n ) stands for the linear interpolation in I n . In this work, only one linear element is used in time within each space-time slab. More details can be found in [START_REF] Johnson | Discontinuous galerkin finite element methods for second order hyperbolic problems[END_REF].

We denote always (u h∆t , v h∆t ) for finite element solutions used in this work (obtained by the time discontinuous space-time Galerkin method with code OOFE, MSSMat, CentraleSupélec): in the chapter 2, they stand for reconstructed time continuous solutions based on the obtained time discrete solutions (u n h∆t , v n h∆t ) ; in the chapters 3 and 4, they stand for time discrete solutions (u n h∆t , v n h∆t ) .

Error estimation methods

Finite element analysis involves always different sources of errors that can compromise the validity of the finite element solutions. One major source of errors is introduced by the spatial-temporal discretization (other sources of errors, such as modelling error, user error, etc. [START_REF] Oden | Estimation of modeling error in computational mechanics[END_REF][START_REF] Utku | Solution errors in finite element analysis[END_REF], exist but our introduction here is mainly related to the discretization error). These errors can be lowered by using smaller mesh sizes or higher order polynomial basis functions, etc., according to some criteria that characterize the accuracy of the finite element solutions. Error estimation is an important tool for quantifying and controlling the errors between the exact solutions and the finite element solutions.

There are mainly two types of error estimates that serve very different purposes:

• a priori error estimate is derived before computing numerical solutions. It tells us the order of convergence for a given finite element method, i.e. how fast the error decreases as the mesh size decreases or the interpolation order increases [START_REF] Babuska | Finite elements: an introduction to the method and error estimation[END_REF][START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. When u is regular enough and u h∆t has some regularities in time, an a priori error estimate usually takes the following form:

∥u -u h∆t ∥ ≤ C(h k + (∆t) l ) (1.26)
where the constant C depends on exact solutions u, k > 0, l > 0, ∥•∥ is some time-space norm. When spatial and temporal approximations are improved, the approximate solutions approach the exact ones and the error goes to 0. However the exact solutions are usually unknown. This method has been already largely applied for the wave equation: Jenkins et al. [START_REF] Jenkins | A priori error estimates for mixed finite element approximations of the acoustic wave equation[END_REF] derived optimal a priori error estimates for mixed finite element displacement formulations of the acoustic wave equation; Deka [START_REF] Deka | A priori L ∞ (L 2 ) error estimates for finite element approximations to the wave equation with interface[END_REF] proposed optimal a priori error estimates for both semidiscrete and fully discrete schemes for the wave equation and it is verified with a numerical example; other contributions include [START_REF] Adjerid | A posteriori finite element error estimation for second-order hyperbolic problems[END_REF][START_REF] Bangerth | Finite element approximation of the acoustic wave equation: Error control and mesh adaptation[END_REF][START_REF] Georgoulis | A posteriori L ∞ (L 2 )-error bounds for finite element approximations to the wave equation[END_REF][START_REF] Hughes | Space-time finite element methods for elastodynamics: formulations and error estimates[END_REF][START_REF] Picasso | Numerical study of an anisotropic error estimator in the L 2 (H 1 ) norm for the finite element discretization of the wave equation[END_REF].

• a posteriori error estimates give us a much better idea of the actual errors in a given finite element computation than a priori estimates. They can be used to perform adaptive mesh refinement [START_REF] Süli | Adaptive finite element approximation of hyperbolic problems[END_REF][START_REF] Tie | Adaptive computation for elastic wave propagation in plate/shell structures under moving loads[END_REF]. More specifically, a posteriori error estimators are used to indicate where the error is particularly high, then the mesh is refined in those locations. A new finite element solution is computed, and the process is repeated until a satisfactory error tolerance is Introduction reached. A posteriori error estimates usually take the form:

∥u -u h∆t ∥ ≤ η = N-1 ∑ n=0 ∑ S∈T h (η n S ) 2 1 2
(1.27)

for dynamic problem. Here η n S (u h∆t ) is a quantity related to u h∆t at t n and in element S ∈ T h . It is called element estimator or indicator. In contrast to a priori estimates, it is used for a practical a posteriori assessment of the accuracy of a computed finite element solution.

One may formulate the following properties describing an optimal a posteriori error estimate [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] Vohralík | A posteriori error estimates for efficiency and error control in numerical simulations[END_REF]:

1. guaranteed upper bound: ensure that eq. (1.27) holds; 2. local efficiency: ensure that the local estimator represents a lower bound for the actual errors, up to a multiplicative constant, i.e. there exist a constant C > 0 such that :

∑ S∈T h (η n S ) 2 1 2 ≤ C ∥u -u h∆t ∥ I n , ∀n ∈ [0, N -1] (1.28)
With these two properties it is known that the error estimator should yield guaranteed and sharp upper and lower bounds of the actual errors.

3. asymptotic accuracy: ensure that the effectivity index:

I eff = η ∥u -u h∆t ∥ (1.29)
i.e. the ratio of the estimated and actual error, goes to one as h and ∆t decrease;

4. robustness: guarantee the three preceding properties independently of the parameters of the problem and of their variations;

5. small evaluation cost: ensure that the computational cost needed for the evaluation of the estimators η n S should be much smaller than the cost required to obtain the approximate solution itself;

A variety of methods are developed for a posteriori error estimates (see [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF] for a review of main methods for error estimation): element residual method [START_REF] Babuška | A posteriori error estimates for the finite element method[END_REF][START_REF] Demkowicz | Toward a universal hp adaptive finite element strategy, part 1. constrained approximation and data structure[END_REF][START_REF] Díez | A posteriori error estimation for standard finite element analysis[END_REF], recovery based methods by Zienkiewicz and Zhu [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineerng analysis[END_REF], the constitutive relation error method by Ladevèze [START_REF] Ladevèze | Constitutive relation error estimators for time-dependent non-linear fe analysis[END_REF][START_REF] Ladeveze | Error estimate procedure in the finite element method and applications[END_REF], goal-oriented dual weighted method [START_REF] Díez | Error estimation and quality control[END_REF][START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF], etc. Most of these methods have been successfully applied for error control of elliptic and parabolic problems [START_REF] Bank | Some a posteriori error estimators for elliptic partial differential equations[END_REF][START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF][START_REF] Bieterman | The finite element method for parabolic equations[END_REF][START_REF] Cottereau | Strict error bounds for linear and nonlinear solid mechanics problems using a patch-based flux-free method[END_REF][START_REF] Cottereau | Fast r-adaptivity for multiple queries of heterogeneous stochastic material fields[END_REF][START_REF] Eriksson | Adaptive finite element methods for parabolic problems i: A linear model problem[END_REF][START_REF] Eriksson | Adaptive finite element methods for parabolic problems iv: Nonlinear problems[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF]. However, there are relatively less studies on the error control of finite element methods for second-order hyperbolic problems [START_REF] Adjerid | A posteriori finite element error estimation for second-order hyperbolic problems[END_REF][START_REF] Bangerth | Finite element approximation of the acoustic wave equation: Error control and mesh adaptation[END_REF][START_REF] Steffens | A simple strategy to assess the error in the numerical wave number of the finite element solution of the helmholtz equation[END_REF][START_REF] Süli | A posteriori error analysis and global error control for adaptive finite element approximations of hyperbolic problems[END_REF].

For example, Georgoulis et al. [START_REF] Georgoulis | A posteriori L ∞ (L 2 )-error bounds for finite element approximations to the wave equation[END_REF] derived a posteriori error bounds of residual type in the L ∞ (L 2 )-norm for for the space-discrete problem and the fully-discrete implicit finite element method of the linear wave equation, using elliptic reconstruction and space-time reconstruction technique. Johnson [START_REF] Johnson | Discontinuous galerkin finite element methods for second order hyperbolic problems[END_REF] proved a priori and a posteriori error estimates for a finite element method for linear second-order hyperbolic equations based on a finite element discretization with discontinuous Galerkin methods. Only the error due to space discretization is considered and the error estimates are derived using elliptic reconstruction. Picasso [START_REF] Picasso | Numerical study of an anisotropic error estimator in the L 2 (H 1 ) norm for the finite element discretization of the wave equation[END_REF] proposed an anisotropic a posteriori error estimate. It was derived for a finite element discretization of the wave equation in two space dimensions and numerical results on adapted meshes indicated that the error estimator slightly underestimates the true error. The works of Aubry et al. [START_REF] Aubry | Adaptive strategy for transient/coupled problems applications to thermoelasticity and elastodynamics[END_REF] were based on a space-time Galerkin formulation for elastodynamics. It provides a variational formulation of the error with respect to the residual. The adjoint state gives an upper bound of the error by the norm of the residuals.

Objectives

The aim of our work is developing tools of error estimation for numerical solutions of wave equation in heterogeneous media.

Firstly, a residual-type a posteriori error representation is proposed for finite element solutions of the elastodynamic equation, based on works of Ibrahima [START_REF] Ibrahima | Estimation d'erreur pour des problèmes de propagation d'ondes en milieux élastiques linéaires hétérogènes[END_REF] for the error estimation of the elastodynamic and acoustic equations. In fact, Vohralík [START_REF] Vohralík | A posteriori error estimates for efficiency and error control in numerical simulations[END_REF] has derived a posteriori error estimates in an energy norm for the first-order (in time) heat equation. Ibrahima attempted to extend his idea to the elastodynamic problem but some difficulties appear when dealing with the second-order time derivative term in time. We continued and prolonged his efforts and developed an explicit error estimation in energy norm for the elastodynamic equation.

Secondly, we plan to quantify HF wave propagation in the weak coupling limit, where strong interactions between wave fronts and the heterogeneous medium generate complex propagation phenomena. This is the primary objective in this work. As discussed in the section 1.2, the energy or more generally quadratic quantities are the most relevant observables in the high frequency limit, and the radiative transfer equation for the energy density arises in the weak coupling regime. We propose that numerical errors can be evaluated with radiative transfer equation in terms of energy quantities of numerical waves. The main objective is defining a new evaluation method of numerical errors of wave fields with this new equation, and determining the distribution of mesh size h related to wave length λ and l c with a more relevant criterion in terms of energy quantities.

In brief, instead of using the classical error method for highly oscillating fields based on the wave equation for u h∆t (fig. 1.6, dashed line), the new errors in terms of energy quantities of wave fields W ε [u h∆t ] (W ε and W denote respectively the Wigner transform and its weak limit, i.e. the Wigner measure) are used to build a posteriori error analysis based on the radiative transfer equation (fig. 1.6, solid line). Difficulties are expected in the aspect of the computation of Wigner transform because of its dependance on four variables (t, x; ω, k) and the complexity of radiative transfer equation for high dimensions.

It is worth noticing that errors quantified by these two methods have different sources. The first method is classical, only the modelling by the elastic wave equation is considered, and we study the discretization error arising in the finite element solutions of this model. to use the radiative transfer model for evaluating the errors of numerical finite element solutions of the wave equation. Thus the defined errors can be considered including both the discretization error and "modelling error", which is the gap between the averaged energy quantities of the exact solutions of the wave equation in a set of "similar" cases and the weak limit defined by the radiative transfer equation for that set.

The manuscript is organized as follows.

In the chapter 2, an explicit a posteriori residual based error upper bound is developed theoretically and numerically for the elastodynamic equation. Basically two main ideas are developed: the secondorder elastodynamic equation is transformed firstly to a first-order hyperbolic system; and the residual method of a posteriori error estimates is exploited with a series of field reconstructions in time and in space. A numerical application with 1D homogeneous and heterogeneous media is given.

Then we develop a residual error estimator for the high frequency wave equation in the weak coupling regime based on radiative transfer modelling in terms of energy quantities. It is realized in homogeneous media in the chapter 3 and heterogeneous media in the chapter 4. This part is the major contribution in this PhD work.

In the chapter 3, we present in detail the Wigner transform and its interest for the analysis of wave equations. We show how it can be used to determine the phase-space energy densities in the high frequency limit and to derive the transport equation in homogeneous media in terms of Wigner measure, the weak limit of the Wigner transform of wave fields. Then a residual type error is based on this equation and it is validated numerically with different mesh sizes in 1D media.

In the chapter 4, following the same idea as the chapter 3, the radiative transfer equation in heterogeneous media is presented. In 1D media, the analytical solutions can be found and the propagation phenomena are introduced and analyzed. A discussion on strong localization is proposed to analyze the numerical results in 1D and to verify the radiative transfer regime. Then the residual errors are defined based on this equation in terms of Wigner transform of approximate solutions of the wave equation. Considering the numerical fluctuations in the calculation of Wigner transform in heterogeneous media, the convolution properties of Wigner transform are used to define a filtered error. Numerical implementation in 1D random media allows to validate these two residual errors. We find that the second one with convolution improves greatly the results of the first one. Some conclusions and perspectives are finally drawn in the chapter 5.

Chapter 2

Explicit a posteriori error estimation for elastodynamic equation This chapter introduces a new a posteriori error estimation for numerical solutions of a secondorder hyperbolic wave equation, based on the work of Ibrahima [START_REF] Ibrahima | Estimation d'erreur pour des problèmes de propagation d'ondes en milieux élastiques linéaires hétérogènes[END_REF] and Vohralík [START_REF] Vohralík | A posteriori error estimates for efficiency and error control in numerical simulations[END_REF]. A key feature is the use of the residual method and the development of a series of reconstructions in time and in space with respect to different regularities required by corresponding ingredients of the obtained error bound.

At first, different reconstructions in time and in space are introduced and discussed in the section 2.1. Then the section 2.2 gives the detailed proof of the obtained estimator. In the section 2.3, numerical results for wave propagation in 1D media are presented for numerical solutions with uniform mesh. Specifically, the behaviour of exact errors and the effectivity of the error bound are studied numerically.

Reconstructions in time and in space of finite element wave solutions

Reconstruction methods have been widely used to construct a posteriori error estimates for finite element approximations. For instance, Mozolevski and Prudhomme [START_REF] Mozolevski | Goal-oriented error estimation based on equilibratedflux reconstruction for finite element approximations of elliptic problems[END_REF] proposed a goal-oriented error estimation in finite element approximations of second-order elliptic problems that combines the dual weighted residual method and equilibrated-flux reconstruction methods for the primal and dual problems. The ZZ-type error estimators due to Zienkiewicz and Zhu [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 1: The recovery technique[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF] are based on reconstruction of an improved stress and define the error as the difference between this stress and the one computed by the standard finite element procedure. The constitutive relation error method proposed by Ladevèze [START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF] consisted in post-processing the finite element solution in order to construct an admissible displacement-stress pair and then defining error bounds by quantifying the non-verification of the constitutive relations.

In this chapter, we consider the case where approximate solutions (u n h∆t , v n h∆t ) at all instant t n are calculated by the finite element method. It means that they are discrete in time (fig. 2.1a). Since time integration properties and some regularities in time are required to deduce the error bound, reconstructions in time can be applied to obtain solutions inside each I n with some regularities. In general, the continuity in space of the exact solution is maintained in the finite element analysis, i.e. u n h∆t and v n h∆t are piecewise polynomial and continuous in space. However, the derivative of finite element solutions in space, and consequently the finite element stress vectors are not continuous across element interfaces [START_REF] Vohralík | A posteriori error estimates for efficiency and error control in numerical simulations[END_REF]. Furthermore, the point-wise strong equilibrium is generally not verified (fig. 2.1b), i.e.

σ σ σ (u n h∆t ) / ∈ H(div, Ω), ρ∂ t v n h∆t -∇ x • σ σ σ (u n h∆t ) ̸ = f (2.1)
where σ σ σ (•) stands for the stress tensor related to a displacement field. Compared to u n h∆t , σ σ σ (u n h∆t ) has the same regularity in time but a reduced regularity in space. Consequently, according to the requirement of regularities by the proposed a posteriori error estimator (section 2.2), some reconstructions in time and in space are proposed here:

• two reconstructions in time for displacement and velocity fields, denoted respectively by u h∆t (t, x) and v h∆t (t, x);

• one reconstruction for stress fields in space, denoted by σ σ σ h∆t (t, x).

These reconstructions are also appropriately chosen in order to simplify or eliminate some terms in the proposed error estimator. We give the detailed definitions of reconstructions that are used in the sections 2.2 and 2.3 as follows.

Firstly, u h∆t and v h∆t satisfy always: ∀n ∈ [0, N -1],

u h∆t (t n , •) = u n h∆t , v h∆t (t n , •) = v n h∆t (2.2)
with initial conditions: u h∆t (0, •) = u 0 and v h∆t (0, •) = v 0 . According to the minimum required regularities for our error estimation, (u h∆t , v h∆t ) should satisfy in time:

u h∆t ∈ C 0 ((0, T )), i.e. C 0 in time (2.3a) v h∆t ∈ C 1 ((0, T )) ∩ ∏ n C 2 (I n ) i.e.
w(t) := (t -t n+1 ) 2 (2t + t n+1 -3t n ) (∆t) 3 w n + (t -t n ) 2 (-2t -t n + 3t n+1 ) (∆t) 3 w n+1 + (t -t n+1 ) 2 (t -t n ) (∆t) 2 ∂ t w n + (t -t n ) 2 (t -t n+1 ) (∆t)
a n h∆t :=            v n+1 h∆t -v n h∆t ∆t for n = 0 v n+1 h∆t -v n-1 h∆t 2∆t for n ∈ [1, N -1] v n h∆t -v n-
) (Ω) + (σ σ σ (u n h∆t ), ε ε ε(w h )) (Ω) = (f n , w h ) (Ω) (2.6)
The computational cost of acceleration in the first v-reconstruction eq. (2.5) is ignorable compared to that in the second one eq. (2.6). Indeed, the computational expense in the latter results mainly from the inverse of the mass matrix. In explicit dynamic analysis, it is at the same cost with the computation of finite element solutions. This is efficient if a diagonalized mass matrix is used. But its influence on reconstructions is not studied here.

Then two definitions are also proposed for the displacement reconstruction u h∆t : h∆t by equilibrium:

Definition 2.
ρ∂ t v n+ 1 2 h∆t , w h ) (Ω) + σ σ σ (u n+ 1 2 h∆t ), ε ε ε(w h ) (Ω) = f n+ 1 2 , w h (Ω)
;

(2.8)

3. reconstruct quadratic u h∆t with the three fields (u n h∆t , u n+1 h∆t , u

n+ 1 2 h∆t ) in each I n .
It is noted that in eq. (2.8), the computation of u h∆t requires to invert the stiffness matrix, which is much more expensive than the inverse of the mass matrix, especially for large models.

Finally, considering the reconstructions in space, since ∇ x • σ σ σ (u h∆t ) is not locally in equilibrium as we introduced before, we propose the following stress reconstruction: Definition 2.5 (Equilibrated σ -reconstruction). Let v h∆t be the reconstructed velocity field (C 1 in time, cubic in each time interval I n ), we call equilibrated σ -reconstruction any function σ σ σ h∆t constructed from v h∆t which satisfies:

σ σ σ h∆t (t, x) ∈ E H(div,Ω),(T ) (ρ∂ t v h∆t -∇ x • σ σ σ h∆t , 1) (S) = (f, 1) (S) ∀t ∈ (0, T ), ∀ S ∈ T h (2.9)
where E H(div,Ω),(T ) denotes space of functions σ σ σ such σ σ σ (t, •) belongs to H(div, Ω) in space and t → σ σ σ (•, x) is in L 2 (0, T ). We define also the reconstructed time derivative of stress:

δ δ δ h∆t (t, x) := ∂ t σ σ σ h∆t , δ δ δ h∆t ∈ E H(div,Ω),(T ) (2.10)
More precisely, σ σ σ h∆t and ∂ t v h∆t have the same regularities in time, i.e. C 0 in time and quadratic in each time interval. If f is not locally quadratic in time intervals, we can take its quadratic interpolation as approximation. Note that eq. (2.9) is only a weak form of the equilibrium equation because it holds only for mean values on each mesh element. Thus σ σ σ h∆t (t, x) is not unique with this definition. We propose the following two reconstructions of σ σ σ h∆t that satisfies the definition 2.5: Definition 2.6 (MRE σ -reconstruction). We call the MRE (Minimum Regularity Ensured) σreconstruction any σ σ σ h∆t , that is continuous and piecewise linear in space and obtained by the following strategies. As ∇ x • σ σ σ h∆t is constant in each S, by eq. (2.9), we have:

∇ x • σ σ σ h∆t = 1 V (S) S (ρ∂ t v h∆t -f) dx (2.11)
where V (S) is the volume of element S. Given the constant derivative of σ σ σ h∆t in each element, one more parameter, such as a Neumann-type boundary condition, is needed so as to determine a continuous and piecewise linear reconstructed stress. σ σ σ h∆t is then determined by minimizing the residual errors.

Definition 2.7 (SA σ -reconstruction). We call SA (Statically Admissible) σ -reconstruction the statically admissible stress construction proposed by Ladevèze and Pelle [START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF], namely it should verify:

∀t, ρ∂ t v h∆t -∇ x • σ σ σ h∆t = f, in S (2.12)
It is known that with the Ladevèze method, at each time step t n , a SA stress field σ σ σ h∆t (t n ) is reconstructed using the weakly equilibrated finite element stress solutions σ σ σ (u n h∆t ). In our case, in order to obtain a stress field σ σ σ h∆t (t) that is continuous in the whole time and quadratic in each time interval, a SA stress field σ σ σ h∆t (t n+ 1 2 ) should be reconstructed and for this the equilibrated ureconstruction at t n+ 1 2 is used. We recall that the latter is based on the equilibrated v-reconstruction according to the definition 2.4. We will see in the section 2.2 that the use of the SA σ -reconstruction allows to cancel one of the indicators and simplify the calculations of the residual error estimators.

The table 2.1 summarizes our propositions of reconstructions presented in this section. They all satisfy the requirement in the obtained error bound and we can make choice according to our demands in numerical applications. In this section, we present the derivation of an a posteriori error upper bound in an energy norm between the unknown exact solutions (u, v) and the reconstructed solutions (u h∆t , v h∆t ). Only homogeneous Dirichlet boundary conditions are considered here for simplicity.

The main result of the explicit a posteriori error upper bound in energy norm is given here: Theorem 2.8. (Error upper bound for elastodynamic problems) Let (u h∆t , v h∆t , σ σ σ h∆t ) a group of reconstruction fields chosen in the table 2.1, we have

∥ v -v h∆t ∥ 2 M(Ω),∞(T ) + ∥ u -u h∆t ∥ 2 K(Ω),∞(T ) ≤ 2∥H R,DF ∥ ∞(T ) + 2H R,DF (0) + 2 √ T ∥H V N ∥ (T ) + 2 √ T ∥H DE,GV ∥ (T ) + √ 2 ∥H R,DF ∥ (T ) ∥H V N ∥ (T ) (2.13) where ∥ v-v h∆t ∥ 2 M(Ω),∞(T ) + ∥ u-u h∆t ∥ 2 K(Ω),∞(T ) = sup t∈[0,T ] ∥ v-v h∆t ∥ 2 M(Ω) (t)+ ∥ u-u h∆t ∥ 2 K(Ω) (t).
The three time dependent residual functions used above are defined below.

The lack of equilibrium function H R,DF (t) is defined as:

H R,DF (t) = ∑ S∈T h (η R,S (t) + η DF,S (t)) 2 1 2
(2.14) with:

η R,S (t) = C p h S C K √ λ min ∥f(t, x) -ρ∂ t v h∆t (t, x) + ∇ x • σ σ σ h∆t (t, x)∥ (S) (2.15) η DF,S (t) = 1 √ λ min ∥σ σ σ h∆t (t, x) -σ σ σ (u h∆t )(t, x)∥ (S) (2.16)
The lack of equilibrium rate function H DE,GV (t) is defined as:

H DE,GV (t) = ∑ S∈T h (η DE,S (t) + η GV,S (t)) 2 1 2
(2.17) with:

η DE,S (t) = C p h S C K √ λ min -∂ t f(t, x) + ρ∂ 2 t v h∆t (t, x) -∇ x • δ δ δ h∆t (t, x) (S) (2.18)
η GV,S (t) = 1 √ λ min ∥-δ δ δ h∆t (t, x) + σ σ σ (∂ t u h∆t (t, x))∥ (S) (2.19)
The velocity residual H V N (t) is defined as:

H V N (t) = ∥∂ t u h∆t (t, x) -v h∆t (t, x) ∥ K(Ω) (2.20)
where h S is the mesh size of element S, λ min is the smallest eigenvalue of the fourth-order tensor C, C p is the constant in the Poincaré's inequality (theorem A.5) and C K is the constant in the Korn inequality (theorem A.6).

In H R,DF , η R,S represents the errors related to equilibrium equation and η DF,S represents the errors between the reconstructed stress and the stress of reconstructed displacement. η R,S vanishes in the case of SA σ -reconstruction according to the definition 2.7. ∥H R,DF ∥ (0) allows involving the initial conditions. H DE,GV results from the residual of time derivative of equilibrium. H V N evaluates the nonconformity between time derivative of the reconstructed displacement and the reconstructed velocity, which equals to 0 only at t n for our choices of reconstructions.

As a preliminary step to the proof of this theorem, we first define a residual operator and prove a lemma related to it. Definition 2.9 (Residual operator R h∆t ). We define the following residual operator R h∆t such that ∀t ∈ (0, T ):

(R h∆t (t), w)

(Ω) = (f(t) -ρ∂ t v h∆t (t), w) (Ω) -(σ σ σ (u h∆t )(t), ε ε ε(w)) (Ω) (2.21)
Explicit a posteriori error estimation for elastodynamic equation where w is considered as an H 1 0 function in space.

According to the regularities of ∂ t v h∆t and u h∆t , R h∆t is C 0 in time and C 1 in each time interval. Its time derivative ∂ t R h∆t (L 2 in time) can be obtained by:

(∂ t R h∆t , w) (Ω) (t) = (∂ t f -ρ∂ 2 t v h∆t , w) (Ω) (t) -(σ σ σ (∂ t u h∆t ), ε ε ε(w)) (Ω) (t) (2.22)
Otherwise, eq. ( 2.21) can be written as:

(R h∆t , w) (Ω) = (f -ρ∂ t v h∆t + ∇ x • σ σ σ (u h∆t ), w) (Ω) + ∑ S∈T h ⟨σ σ σ (u h∆t ) • n, w⟩ (∂ S) (2.23)
where the last term gives rise to the stress vector jumps at the interfaces of elements. Thus this operator involves both the volume and surface residuals. Also, by adding and substracting the equilibrated stress reconstruction σ σ σ h∆t into the definition of residual R h∆t in eq. ( 2.21), we have

(R h∆t , w) (Ω) = (f -ρ∂ t v h∆t + ∇ x • σ σ σ h∆t , w) (Ω) + (σ σ σ h∆t -σ σ σ (u h∆t ), ε ε ε(w)) (Ω) (2.24)
It can be seen that this operator contains also the residual of different reconstructions in time.

Finally, according to the weak formulation, we have:

(R h∆t , w) (Ω) (t) = (ρ∂ t (v -v h∆t ), w) (Ω) (t) -(σ σ σ (u -u h∆t ), ε ε ε(w)) (Ω) (t) (2.25)
The definition of ∂ t R h∆t in eq. (2.22) and the property in eq. (2.25) will be used in the derivation of the theorem 2.8. Lemma 2.10 (Upper bound of residual R h∆t ). ∀t ∈ (0, T ), we have:

(R h∆t , w) (Ω) (t) ≤ H R,DF (t) ∥w∥ K(Ω) (t) (2.26) and t 0 (R h∆t , ∂ t w) (Ω) dτ ≤ ∥H DE,GV ∥ (t) ∥w∥ K(Ω),(t) + H R,DF (t) ∥w∥ K(Ω) (t) + H R,DF (0) ∥w∥ K(Ω) (0)
(2.27) where w is H 1 0 in space, C 0 in time and C 1 in each time interval.

Proof. We prove firstly the equation (2.26) based on eq. (2.24). According to the definition 2.5 of σ σ σ h∆t and assuming that w S is the average of w in element S, we have:

(ρ∂ t v h∆t -∇ x • σ σ σ h∆t , w S ) (S) = (f, w S ) (S) , ∀ S (2.28)
Decomposing eq. (2.24) in space into each element S and introducing eq. (2.28) into eq. ( 2.24), we can deduce that:

(R h∆t , w) (Ω) = ∑ S∈T h (f -ρ∂ t v h∆t + ∇ x • σ σ σ h∆t , w -w S ) (S) + ∑ S∈T h (σ σ σ h∆t -σ σ σ (u h∆t ), ε ε ε(w)) (S) (2.29)
For the first term in eq. (2.29), using respectively Cauchy-Schwarz inequality (theorem A.7), Poincaré's inequality, Korn's inequality and that λ min is the smallest eigenvalue of the fourth-order tensor C:

(f -ρ∂ t v h∆t + ∇ x • σ σ σ h∆t , w -w S ) (S) ≤ ∥f -ρ∂ t v h∆t + ∇ x • σ σ σ h∆t ∥ (S) ∥w -w S ∥ (S) ≤ C p h S ∥f -ρ∂ t v h∆t + ∇ x • σ σ σ h∆t ∥ (S) ∥w∥ 1(S) ≤ C p h S C K ∥f -ρ∂ t v h∆t + ∇ x • σ σ σ h∆t ∥ (S) ∥ε ε ε(w)∥ (S) ≤ C p h S C K √ λ min ∥f -ρ∂ t v h∆t + ∇ x • σ σ σ h∆t ∥ (S) ∥w∥ K(S) = η R,S ∥w∥ K(S) (2.30)
where η R,S is defined in eq. (2.15). For the second term in eq. (2.29), using again the Cauchy-Schwarz inequality and the definition of λ min results in:

(σ σ σ h∆t -σ σ σ (u h∆t ), ε ε ε(w)) (S) ≤ ∥σ σ σ h∆t -σ σ σ (u h∆t )∥ (S) ∥ε ε ε(w)∥ (S) ≤ 1 √ λ min ∥σ σ σ h∆t -σ σ σ (u h∆t )∥ (S) ∥w∥ K(S) = η DF,S ∥w∥ K(S) (2.31)
where η DF,S is defined in eq. (2.16). Introducing eq. (2.30) and eq. (2.31) into eq. (2.29), then applying again the Cauchy-Schwarz inequality, we obtain eq. (2.26):

(R h∆t , w) Ω (t) ≤ ∑ S∈T h (η R,S (t) + η DF,S (t)) ∥w∥ K(S) (t) ≤ ∑ S∈T h (η R,S (t) + η DF,S (t)) 2 1 2 ∥w∥ K(Ω) (t) = H R,DF (t) ∥w∥ K(Ω) (t) (2.32)
Now the same idea of proof can be applied for the second result eq. (2.27) of lemma 2.10. Replacing w by ∂ t w, applying integration by parts in time and using eq. (2.22), we have:

t 0 (R h∆t , ∂ t w) (Ω) dτ = - t 0 (∂ t R h∆t , w) (Ω) dτ + (R h∆t , w) (Ω) (t) -(R h∆t , w) (Ω) (0) (2.33)
Explicit a posteriori error estimation for elastodynamic equation Adding and substracting the reconstructed time derivative of stress (eq. (2.10)) into the first term of eq. (2.33) results in:

- t 0 (∂ t R h∆t , w) (Ω) dτ = t 0 (-∂ t f + ρ∂ 2 t v h∆t -∇ x • δ δ δ h∆t , w) (Ω) dτ - t 0 (δ δ δ h∆t -σ σ σ (∂ t u h∆t ), ε ε ε(w)) (Ω) dτ
(2.34) It can be observed that we can proceed with the same steps of proof with eq. (2.24). We have:

-(∂ t R h∆t , w) (Ω) (t) ≤ H DE,GV (t) ∥w∥ K(Ω) (t) (2.35)
where H DE,GV is defined in eq. (2.17). Finally, introducing eq. ( 2.35) into eq. (2.33), using the Cauchy-Schwartz inequality to obtain time-space norm, then applying eq. ( 2.26) for the last two terms of eq. ( 2.33), we deduce:

t 0 (R h∆t , ∂ t w) (Ω) dτ ≤ ∥H DE,GV ∥ (t) ∥w∥ K(Ω),(t) + H R,DF (t) ∥w∥ K(Ω) (t) + H R,DF (0) ∥w∥ K(Ω) (0) (2.36)
Now we start to prove theorem 2.8:

Proof. (Theorem 2.8)
∀t, we have:

1 2 ∥ v -v h∆t ∥ 2 M(Ω) (t) + 1 2 ∥ u -u h∆t ∥ 2 K(Ω) (t) = 1 2 t 0 d dt ∥v -v h∆t ∥ 2 M(Ω) + d dt ∥u -u h∆t ∥ 2 K(Ω) dτ + 1 2 ∥v -v h∆t ∥ 2 M(Ω) (0) + 1 2 ∥u -u h∆t ∥ 2 K(Ω) (0) = t 0 (ρ∂ t (v -v h∆t ), v -v h∆t ) (Ω) dτ + t 0 (σ σ σ (u -u h∆t ), ε ε ε(∂ t (u -u h∆t ))) (Ω) dτ
(2.37) The last two terms related to the initial conditions vanish according to hypothesis for reconstruction fields at t = 0: u h∆t (0, •) = u 0 and v h∆t (0, •) = v 0 . Then by adding and subtracting ∂ t u h∆t in the first term of the last line of eq. (2.37) and using the definition 2.9 of R h∆t with w = ∂ t (u -u h∆t ), we get:

1 2 ∥ v -v h∆t ∥ 2 M(Ω) (t) + 1 2 ∥ u -u h∆t ∥ 2 K(Ω) (t) = t 0 (ρ∂ t (v -v h∆t ), ∂ t u h∆t -v h∆t ) (Ω) dτ + t 0 (ρ∂ t (v -v h∆t ), v -∂ t u h∆t ) (Ω) dτ + t 0 (σ σ σ (u -u h∆t ), ε ε ε(∂ t (u -u h∆t ))) (Ω) dτ = t 0 (ρ∂ t (v -v h∆t ), ∂ t u h∆t -v h∆t ) (Ω) dτ :=(a) + t 0 (R h∆t , ∂ t (u -u h∆t )) (Ω) dτ :=(b) (2.38)
The term (b) is obtained according to eq. (2.25). Then the two terms (a) and (b) are treated separately in the following.

Firstly, for the term (a), using eq. (2.25) by taking w = (∂ t u h∆t -v h∆t ), we have:

(a) = t 0 (R h∆t , ∂ t u h∆t -v h∆t ) (Ω) dτ + t 0 (σ σ σ (u h∆t -u), ε ε ε(∂ t u h∆t -v h∆t )) (Ω) dτ (2.39)
Now the first equation (2.26) of lemma 2.10 can be applied for the first term of eq. (2.39):

t 0 (R h∆t , ∂ t u h∆t -v h∆t ) (Ω) dτ ≤ t 0 H R,DF (τ) ∥∂ t u h∆t -v h∆t ∥ K(Ω) (τ)dτ ≤ ∥H R,DF ∥ (t) ∥∂ t u h∆t -v h∆t ∥ K(Ω),(t) ≤ ∥H R,DF ∥ (T ) ∥H V N ∥ (T ) (2.40)
where H V N is defined in eq. (2.20). The second term of eq. (2.39) is treated firstly by Cauchy-Schwartz inequality:

t 0 (σ σ σ (u h∆t -u), ε ε ε(∂ t u h∆t -v h∆t )) (Ω) dτ ≤ ∥u -u h∆t ∥ K(Ω),(t) ∥∂ t u h∆t -v h∆t ∥ K(Ω),(t) = ∥u -u h∆t ∥ K(Ω),(t) ∥H V N ∥ (t) (2.41)
Using the definition of L ∞ norm, it can be derived that:

∥u∥ 2 K(Ω),(t) = t 0 ∥u∥ 2 K(Ω) (τ)dτ ≤ t 0 sup τ ∥u∥ 2 K(Ω) (τ) dτ = t sup τ∈[0,t] ∥u∥ 2 K(Ω) (τ) = t ∥u∥ 2 K(Ω),∞ (t)
(2.42) Apply this result for ∥u -u h∆t ∥ K(Ω),(t) in eq. (2.41):

t 0 (σ σ σ (u h∆t -u), ε ε ε(∂ t u h∆t -v h∆t )) (Ω) ≤ √ t ∥u -u h∆t ∥ K(Ω),∞(t) ∥H V N ∥ (t) ≤ √ T ∥H V N ∥ (T ) ∥u -u h∆t ∥ K(Ω),∞(t) (2.43)
Combing eq. (2.40) and eq. (2.43), the temporary upper bound is found for (a):

(a) ≤ ∥H R,DF ∥ (T ) ∥H V N ∥ (T ) + √ T ∥H V N ∥ (T ) ∥u -u h∆t ∥ K(Ω),∞(T ) (2.44)
Secondly, for the term (b), taking w = uu h∆t in the second equation (2.27) of lemma 2.10, we have:

(b) ≤ ∥H DE,GV ∥ (t) ∥u -u h∆t ∥ K(Ω),(t) + H R,DF (t) ∥u -u h∆t ∥ K(Ω) (t) + H R,DF (0) ∥u -u h∆t ∥ K(Ω) (0) ≤ √ T ∥H DE,GV ∥ (T ) + ∥H R,DF ∥ ∞(T ) + H R,DF (0) ∥u -u h∆t ∥ K(Ω),∞(T ) (2.45)
Explicit a posteriori error estimation for elastodynamic equation

We have successfully found an upper bound for (a) and (b) in eq. (2.44) and eq. (2.45), then inserting these two inequalities into eq. (2.38), ∀t ∈ (0, T ], we have:

∥ v -v h∆t ∥ 2 M(Ω) (t)+ ∥ u -u h∆t ∥ 2 K(Ω) (t) ≤ A ∥u -u h∆t ∥ K(Ω),∞(t) + B 2 ≤ A ∥v -v h∆t ∥ 2 M,∞(t) + ∥u -u h∆t ∥ 2 K(Ω),∞(t) + B 2 (2.

46)

where:

A = 2∥H R,DF ∥ ∞(T ) + 2H R,DF (0) + 2 √ T (∥H DE,GV ∥ (T ) + ∥H V N ∥ (T ) ) B 2 = 2∥H R,DF ∥ (T ) ∥H V N ∥ (T ) (2.47) As x 2 ≤ ax + b 2 implies x ≤ a + b whenever a, b positive 1
, we find the a posteriori error bound:

∥v -v h∆t ∥ 2 M(Ω),∞(T ) + ∥u -u h∆t ∥ 2 K(Ω),∞(T ) ≤ A + B

Numerical application in 1D homogeneous and heterogeneous media

In this section, we apply numerically the theorem 2.8. Wave propagation in a one-dimensional elastic bar with homogeneous media and two different heterogeneous media is considered here. In these cases, the analytical solutions or reference solutions of the wave equation can be found. We compare numerically the "exact errors" (i.e. errors between analytical or reference solutions and numerical solutions) and the estimated errors defined in eq. (2.13) for numerical solutions with different mesh sizes. The efficiency of the error upper bound is analyzed.

Definition of numerical model and parameters

We consider here a 1D bar with length L = 3.2 m, whose left end is subjected to the Neumann boundary condition with a truncated ricker signal F(t) with a magnitude of 1 kN/m2 and the duration T R of 0.22 ms (fig. 2.3). The right end of the bar is subjected to the free boundary condition. The total analysis time period is [0, T f ] with T f = 1 ms which allows the wave front to propagate a round trip in this bar. 2) a heterogeneous bar with one inclusion, with E = E 2 = E 1 /70 in the inclusion and E = E 1 otherwise, where the length of the inclusion is l = 0.2 m, and its centre is located at L/2 = 1.6 m;

3) a heterogeneous bar with five inclusions, with E = E 2 in the inclusions and E = E 1 otherwise, where the length of each inclusion is l = 0.2 m, and their centres are located respectively at 1, 1.4, 2, 2.4, 2.8 m. For this bar, the observation period is changed to 3T f to ensure a round trip because the wave velocity in inclusions is smaller.

The numerical parameters are given in the tables 2.2, 2.3 and 2.4. In the table 2.2, f c is the cutoff frequency of the truncated ricker signal. The table 2.3 presents the different mesh sizes h and time steps ∆t used for the homogeneous case and the table 2.4 presents those used for both heterogeneous cases. It can be observed that in heterogeneous cases we have more elements per wavelength in the Explicit a posteriori error estimation for elastodynamic equation slow area with E 1 than that with E 2 . Thus the obtained numerical solutions describe waves more precisely in these areas. media load (truncated ricker)

ρ (kg/m 3 ) E (GPa) c (m/s ) T R (ms) f c (kHz) homogeneous 2500 E 1 = 70 c 1 = 5300 0.22 9.5 heterogeneous 2500 E 1 = 70 c 1 = 5300 0.22 9.5 E 2 = 1 c 2 = 632

Analytical solutions of wave equation

In this section, analytical or reference solutions are given for the homogeneous case and the two heterogeneous cases. It allows computing their differences with numerical solutions.

For the homogeneous case, analytical solutions exist:

u(t, x) =              0 for 0 ≤ t ≤ x c 1 c 1 E t-x c 1 0 F(τ)dτ for x c 1 ≤ t ≤ L c 1 c 1 E t-x c 1 0 F(τ)dτ -c 1 E t-L c 1 + x c 1 -L c 1 0 F(τ)dτ for L c 1 ≤ t ≤ 2L c 1 (2.48)
For the heterogeneous bar with one inclusion, analytical solutions for the wave equation can be also found. In fig. 2.5, the reflection and transmission of waves at an interface x = a is given. Since the incident stress waves is defined by σ i = F(tx/c 1 ), using the continuity of the stress at each interface, the stresses for reflected and transmitted waves can be expressed as [START_REF] Achenbach | Wave propagation in elastic solids[END_REF]:

σ r = C r F(t + (x -2a)/c 1 ), σ t = C t F(t -a/c 1 -(x -a)/c 2 ) (2.49)
with C r and C t are respectively the reflection coefficient and the transmission coefficient given by:

C r = c 2 /c 1 -1 c 2 /c 1 + 1 , C t = 2c 2 /c 1 c 2 /c 1 + 1 (2.50)
Then total waves can be calculated with the superposition of reflected and transmitted waves.

For the heterogeneous bar with five inclusions, considering the complexity of analytical solutions under a large number of reflections and transmissions, numerical solutions with a relatively small mesh size h = 1 mm are used as the reference solution.

F (t -x/c 1 ) C t F (t -a/c 1 -(x-a)/c 2 ) (ρ, c 1 ) C r F (t +(x-2a)/c 1 ) x=a (ρ, c 2 )
Fig. 2.5 Incident, reflected and transmitted waves

Numerical exact errors for 1D elastic wave propagation

In this section, we observe firstly the errors between analytical solutions found in the section 2.3.2 and numerical finite element solutions in energy norm in 1D media. Firstly, fig. 2.6 presents an example of numerically calculated velocity fields and local exact errors in energy norm in time-space between approximate solutions and analytical or reference solutions:

e M+K (t, x) = 1 2 ρe 2 v + E(∂ x e u ) 2 (t, x), with e v = v -v h∆t , e u = u -u h∆t (2.51)
It is observed that overall, the errors increase with time and with number of interaction with the interfaces or boundaries. Then the evolution of errors with time:

∥e v ∥ M(Ω) (t) = L 0 1 2 ρe 2 v (t, x)dx 1 2 , ∥e u ∥ K(Ω) (t) = L 0 1 2 E(∂ x e u (t, x)) 2 dx 1 2 (2.52)
is studied so as to observe the error behaviour on the interfaces and the boundary of domain.

In the homogeneous case, ∥e v ∥ M(Ω) (t) and ∥e u ∥ K(Ω) (t) are presented in fig. 2.7b. Firstly, errors are increasing with time of simulation. A strong oscillation observed in [t 1 ,t 2 ] is related to when waves interact with the right boundary (fig. 2.7a). A conversion between errors in two energy norms is obtained, resulting from the conversion of potential energy and kinetic energy. Now the errors for the first heterogeneous case are shown in fig. 2.8b. We observe always the increase of errors with time and strong oscillations during [t 1 ,t 2 ] and [t 3 , T f ]. From fig. 2.8a it can be known that these two time intervals are corresponding to the period when waves interact either with the two boundaries of bar or the interfaces between two materials. During [t 2 ,t 3 ], wave front propagates far from the interfaces or boundaries. Note that during [0,t 1 ], waves are reflected and refracted by the first interface around t = 3 ms but this effect on errors is relative small compared to that in latter times so we neglect here (inset of fig. T (s)

×10 -3 t 1 t 2 t 3 ∥e v ∥ M ∥e u ∥ K (b)
∥e v ∥ M (t) and ∥e u ∥ K(Ω) (t)
Finally, in the highly heterogeneous case with 5 inclusions, it is difficult to identify the instants as what we have done above since waves are largely reflected and refracted here. However, we can still find that oscillation of errors start at t 1 when waves meet the first interface (fig.

T (s) ×10 -3 t 1 ∥e v ∥ M ∥e u ∥ K t 1 (b)
∥e v ∥ M (t) and ∥e u ∥ K(Ω) (t)
Now we focus on our L ∞ -time exact errors in energy norm defined in the error bound:

∥v -v h∆t ∥ 2 M(Ω),∞(T f ) + ∥u -u h∆t ∥ 2 K(Ω),∞(T f ) = ∥e v ∥ 2 M(Ω),∞(T f ) + ∥e u ∥ 2 K(Ω),∞(T f ) (2.53)
Equation (2.53) is computed for all the mesh sizes in homogeneous and heterogeneous media and we present all the results in terms of degree of freedom in fig. 2.10.

It can be observed in fig. 2.10 that in all cases, exact errors are decreasing with smaller mesh size or larger number of degrees of freedom and more precisely. They are inversely proportional to dof. Errors increase with the number of degrees of heterogeneity for solutions with the same degree of freedom.

v ∥ 2 M(Ω),∞(T f ) + ∥e u ∥ 2 K(Ω),∞(T f ) (black line is used as reference with slope -1)

Numerical estimated errors for 1D elastic wave propagation

Now we compute numerically the error bound in eq. (2.13) and verify its accuracy. Four combinations of the reconstructions in time and in space listed in the table 2.1 of the section 2.1 are studied:

• case1: basic u-reconstruction + basic v-reconstruction + MRE σ -reconstruction; • case2: basic u-reconstruction + equilibrated v-reconstruction + MRE σ -reconstruction; • case3: equilibrated u-reconstruction + equilibrated v-reconstruction + MRE σ -reconstruction; • case4: equilibrated u-reconstruction + equilibrated v-reconstruction + SA σ -reconstruction.
Recall that the SA σ -reconstruction requires both the equilibrated u-reconstruction and the equilibrated v-reconstruction. In the rest of four combinations, we choose three of them here (case1-3).

The error bound is firstly computed in homogeneous media.

In fig. 2.11 we compared numerical exact errors obtained in the section 2.3.3 and estimated errors with the four cases presented above.

Firstly, it is illustrated in fig. 2.11a that estimated errors with four different reconstructions are decreasing with dof, as we expected. They are close to each other. The second case gives the smallest value of errors. Its computational cost is relatively small, especially compared to the last two cases. Secondly, compared to the numerical exact errors, the obtained error bounds are much larger. On the one hand, it can be understood by the existence of a significant constant coefficient before the estimated errors; on the other hand, their rate of convergence is different. In effect, in fig. 2.11b we find that the effectivity index of error bound (i.e. ratio of the estimated and exact errors, see eq. (1.29)) increases with the computational effort. Thus this estimator has a weak asymptotic accuracy.

√ T H V N (T ) ) case1 case2 case3 case4 exact (b)
Fig. 2.12 A log-log plot of two terms of error bound: (a) 2

√ T ∥H DE,GV ∥ (T ) (b) 2 √ T ∥H V N ∥ (T )
Now we analyze the major terms that contribute to total estimated errors. In fig. 2.12a, we find 2 √ T ∥H DE,GV ∥ (T ) is dominant in estimated errors. Recall that:

η DE,S = C p h S C K √ λ min -∂ t f + ρ∂ 2 t v h∆t -∇ x • δ h∆t (S) , η GV,S = 1 √ λ min ∥-δ h∆t + σ (∂ t u h∆t )∥ (S) (2.54)
so it depends on both reconstructions in time and in space. We find that in fact this term contributes most to the total errors due to the nonconformity of the time derivative of equilibrium of reconstruction fields.

In fig. 2.12b, ∥H DE,GV ∥ (T ) also contributes a large part of errors in the cases 3 and 4. Recall that ∥H DE,GV ∥ (T ) = ∥∂ t u h∆tv h∆t ∥ K(Ω),(T ) so it depends mainly on the reconstructions of u and v in time; it is found in fig. 2.12b that the cases 1 and 2 with basic reconstruction for u, v have a better conformity between ∂ t u h∆t and v h∆t inside I n than that in the cases 3 and 4 with equilibrated u-reconstruction and v-reconstruction. Besides, for the cases 3 and 4 same reconstructions (u h∆t , v h∆t ) are used so they give the same results.

Conclusion

In this chapter we derived explicit a posteriori error estimates in a non-natural L ∞ norm for elastic wave propagation in heterogeneous media. The numerical results in the homogeneous case indicate that this estimation gives a fully computable upper bound of exact errors. However, its asymptotic exactness remains to be improved. In the work of Vohralík [START_REF] Vohralík | A posteriori error estimates for efficiency and error control in numerical simulations[END_REF], similar reconstructions were also defined and used for an error bound of the first order (in time) heat equation; a nearly constant effectivity index was obtained in its numerical applications. In fact, for the second-order (in time) hyperbolic equation, our error bound includes not only the residual errors related to the equilibrium (as in [START_REF] Vohralík | A posteriori error estimates for efficiency and error control in numerical simulations[END_REF]) but also the residual errors related to the time derivative equilibrium of reconstructed fields. The latter (H DE,GV ) dominates in the behaviour of the estimated errors. In a future work, more efforts should focus on reconstructions of fields that respect the equilibrium better. As we introduced in the section 1.2.4, the radiative transfer equation arises for high frequency wave propagation in heterogeneous media in the weak coupling limit. It is derived by a two-scale asymptotic expansion of the wave equation in terms of the spatio-temporal Wigner transforms of wave fields. Its particular form in homogeneous media, i.e. the transport equation, is studied in this chapter. It describes the transport of wave energy. Considering the highly oscillating properties of HF wave fields, the classical error estimator that was derived in chapter 2 is no more efficient. The errors based Energy-based error estimation for HF wave in homogeneous media on more relevant quantities, i.e. the energy densities of numerical solutions of wave equation, are studied here with the transport equation.

In the section 3.1, the Wigner transform is introduced in detail. We start with the temporal Wigner transform of a time signal. Its important mathematical properties, especially the interference property and the energy property, are presented and analyzed. Then the scaled spatio-temporal Wigner transform in the high frequency limit is defined and it satisfies also all the properties above. Finally, the discrete Wigner transform is studied and some numerical examples are given.

In the section 3.2, the principal steps for the derivation of the transport equation in terms of Wigner transforms are given. The transition from wave fields to energy densities, from the wave equation model to the energy transport model is presented in detail. Then the theoretical residual errors of finite element solutions of the wave equation based on the transport equation are defined and developed.

In the section 3.3, the numerical application of the defined residual error in a 1D medium is given. Residual errors for finite element solutions with different discretizations in time and in space are compared, that allows validating the obtained error estimator.

Wigner transform applied to the elastic wave equation

In this section, we give a detailed introduction of the Wigner transform. In the high frequency regime, it is used to derive the transport equation in homogeneous media and the radiative transfer equation in heterogeneous media (chapter 4). Some examples and numerical illustrations are given to clarify its properties and physical meanings.

The Wigner transform is widely used in signal processing [START_REF] Boashash | Time-frequency signal analysis and processing: a comprehensive reference[END_REF][START_REF] Cohen | Time-frequency distributions-a review[END_REF][START_REF] Sandsten | Time-frequency analysis of time-varying signals and non-stationary processes[END_REF] as a time-frequency representation since it provides a better simultaneous time-frequency resolution compared to the other representations, such as the Fourier transform (FT) [START_REF] Bloomfield | Fourier analysis of time series: an introduction[END_REF][START_REF] Bracewell | The Fourier transform and its applications[END_REF], the short-time Fourier transform (STFT) [START_REF] Auger | Time-frequency toolbox[END_REF][START_REF] Baba | Time-frequency analysis using short time fourier transform[END_REF], etc.

For time-space dependent functions, the definition of spatio-temporal Wigner transform is given by: Definition 3.1 (Spatio-temporal Wigner transform). The spatio-temporal Wigner transform of two time-space dependent vector functions F(t, x) and G(t, x), (t, x) ∈ R × R d , is a second-order tensor defined by:

W[F, G](t, x; ω, k) = 1 (2π) 1+d R×R d F t - τ 2 , x - y 2 ⊗ G * t + τ 2 , x + y 2 e i(τω+y•k) dτdy (3.1)
where d is the spatial dimension and the superscript * denotes the complex conjugation operator.

Temporal Wigner transform in signal processing

The spatio-temporal Wigner transform (3.1) depends on four variables, two in the time-space domain and two in the (ω, k)-phase space. In order to simplify the introduction, we give the definition of temporal Wigner transform for two scalar time signals and introduce its properties in the following. The spatio-temporal Wigner transform satisfies also these properties. Definition 3.2 (Temporal Wigner transform). The (cross) Wigner transform of two time signals F(t) and G(t) is defined as the Fourier transform of their instantaneous correlation function:

W[F, G](t; ω) = 1 2π R F t - τ 2 G * t + τ 2 e iτω dτ (3.2)
The (auto-) Wigner transform of F(t) is denoted by: W

[F](t; ω) := W[F, F](t; ω).
Basically, the temporal Wigner transform tells us how the spectral density changes in time. It has several important properties [START_REF] Sandsten | Time-frequency analysis of time-varying signals and non-stationary processes[END_REF]:

Property 1 (Symmetry property). If F(t) is real-valued, W[F] is an even function of ω. W[F](t; ω) = W[F](t, -ω) (3.
|W[F 1 + F 2 ](t; ω)| 2 = R×R W[F 1 ] t - τ 2 ; ω - υ 2 W[F 2 ] t + τ 2 ; ω + υ 2 dτdυ (3.5)
Property 4 (Integration properties). The cross Wigner transform can be seen as a time distribution of the product of spectrums F Ĝ * or a phase-space description of the product of time signals

FG * . R W[F, G](t; ω)dt = 1 2π F(ω) Ĝ * (ω), R W[F, G](t; ω)dω = F(t)G * (t) (3.6)
Property 5 (Energy density properties). The marginal distribution obtained by integrating over time equals the energy spectrum, while the one obtained by integrating over frequency equals the instantaneous energy.

R W[F](t; ω)dt = 1 2π F(ω) 2 , R W[F](t; ω)dω = |F(t)| 2 (3.7)
The integration of Wigner transform over the time-frequency domain is equal to the total energy of the signal.

R×R W[F](t; ω)dtdω = R |F(t)| 2 dt = 1 2π R F(ω) 2 dω (3.8)
The two equations above mean that Wigner transform can be interpreted as a distribution of energy density in time-frequency domain, which is, however, not always positive for all (t; ω). Property 6 (Convolution property). The Wigner transform of the product of two signals equals the convolution of the Wigner transform of each signal with respect to frequency. The Wigner transform of the convolution of two signals equals the convolution of the Wigner transform of each signal with respect to time.

W[Fh](t; ω) = R W[F](t; υ)W[h](t; ω -υ)dυ (3.9a) W[F * h](t; ω) = R W[F](τ; ω)W[h](t -τ; ω)dτ (3.9b)
where " * " denotes the convolution operator of two time functions.

Property 7 (Interference property). The Wigner transform of a sum of two signals has the following form:

W[F 1 + F 2 ](t; ω) = W[F 1 ](t; ω) + W[F 2 ](t; ω) + 2Re {W[F 1 , F 2 ](t; ω)} (3.10)
where "Re" denotes the real part of a complex-valued quantity. This is due to the fact that Wigner transform is bilinear function according to its definition 3.2. W[F 1 ] and W[F 2 ] are called "auto-terms" and 2Re {W[F 1 , F 2 ](t; ω)} is called "cross-term" or interference term.

Discussion and analysis of the interference property of Wigner transform

The existence of cross-terms of Wigner transform has a great influence on the time-frequency analysis. In this section we study in detail this property by some theoretical and numerical examples.

The following example is given by [START_REF] Hlawatsch | Interference terms in the wigner distribution[END_REF][START_REF] Janssen | Application of the Wigner distribution to harmonic analysis of generalized stochastic processes[END_REF]. We consider a time signal F(t) = F 1 (t) + F 2 (t) defined as follows:

F 1 (t) = A 1 h(t -t 1 )e iω 1 t , F 2 (t) = A 2 h(t -t 2 )e iω 2 t (3.11)
where A 1 and A 2 are real constant amplitudes, h is thought to be a window function concentrated at 0 in time. Thus F 1 and F 2 are respectively concentrated around (t 1 ; ω 1 ) and (t 2 ; ω 2 ) in the (t; ω)-plane.

Using the interference property eq. (3.10), W[F 1 + F 2 ] can be computed analytically. The two auto-terms can be obtained by the convolution property eq. (3.9) and by taking into account the fact that W[e iω 1 t ] = δ (ωω 1 ):

W[F 1 ](t; ω) = A 2 1 W[h](t -t 1 ; ω -ω 1 ), W[F 2 ](t; ω) = A 2 2 W[h](t -t 2 ; ω -ω 2 ) (3.12)
Thus as a time-frequency representation of F 1 and F 2 , W[F 1 ] and W[F 2 ] lie respectively in the surroundings of points (t 1 , ω 1 ) and (t 2 , ω 2 ) in the (t; ω)-plane. For the cross term 2Re {W[F 1 , F 2 ](t; ω)}, by the definition of the cross Wigner transform eq. (3.2) and a time shift of h, we obtain:

2Re {W[F 1 , F 2 ](t; ω)} = 2Re A 1 A 2 1 2π h t -t 1 - τ 2 e iω 1( t-τ 2 ) h * t -t 2 + τ 2 e -iω 2( t+ τ 2 ) e iτω dτ = 2A 1 A 2 cos (ω d t -(ω -ω m )t d ) W[h](t -t m ; ω -ω m ) (3.
13) where the time lag and the frequency lag are denoted by (t d ; ω d ) = (t 2t 1 ; ω 2ω 1 ), the centre time and the centre frequency are denoted by (t m ; ω m ) = ( t 1 +t 2 2 ; ω 1 +ω 2 2 ). Several properties of cross term are observed with this expression. Firstly, the cross term is a modulated version of the original signal F shifted to the midpoint (t m ; ω m ) of the segment connecting (t 1 ; ω 1 ) and (t 2 ; ω 2 ) in the (t; ω)-plane. This result can be generalized to a function composed of N time-frequency component signal, which are localized around N different points in the (t; ω)-plane. In this case, there exist N(N-1) 2 cross terms between each two components. Secondly, the "cosine" function in cross term can be written as:

cos (ω d t -(ω -ω m )t d ) = cos   ω 2 d + t 2 d   ω d ω 2 d + t 2 d t -ω t d ω 2 d + t 2 d   + ω m t d   (3.14)
so it results in oscillations in the direction (ω d ; -t d ) that is orthogonal to the line that connects the auto-terms in the direction (t d ; ω d ); the frequency of this oscillation is simply the distance

t 2 d + ω 2 d
Energy-based error estimation for HF wave in homogeneous media in the (t; ω)-plane between the signal components. Then the peak value of the cross term is twice as large as the product of that of the auto-terms. Finally, we will see with the following example that the "cosine" function is in fact the source of negative values of Wigner transform.

Numerical illustration of Wigner transform of a sum of two gaussian complex-valued signals (fig. 3.1a) is given in the following so as to observe the properties of cross terms explained above. We define h(t) and the parameters of the function F in eq. (3.11) by:

h(t) = e -2t 2 , A 1 = ( 4
π ) 1 4 , (t 1 ; ω 1 ) = (-2; 20); A 2 = ( 8
π ) 1 4 , (t 2 ; ω 2 ) = (2; 40) (3.15)
Figure 3.1b gives the Wigner transform of F(t). It is found that the two auto-terms located at (t 1 ; ω 1 ), (t 2 ; ω 2 ) and the cross term located at their centre (t m ; ω m ) = (0; 30) agree with the theoretical analysis above. We also observe that oscillations of the cross term (along the solid arrowed line) are orthogonal to connecting line of the auto-terms (along the dashed arrowed line). The auto-term W [F 2 ] and the cross term have the same amplitude and quadruple that of the auto-term As expected, it is observed that the cross term contains some negative values, resulting from "cosine" in eq. (3.13). Although the Wigner transform satisfies always energy properties (eq. (3.7) and (3.8)), we cannot interpret it as an exact time-varying energy density estimate. It may be more reasonable to integrate Wigner transform in a period of time or a band of frequency so as to ensure its positivity as a representation of energy (energy density property eq. (3.7)). The integral of autoor cross terms over frequency is called its energy content. Generally, both auto-and cross terms contribute to signal energy [START_REF] Hlawatsch | Interference terms in the wigner distribution[END_REF][START_REF] Qian | Wigner distribution decomposition and cross-terms deleted representation[END_REF]. The oscillations of cross terms suggest that their energy content is small compared to the energy content of auto-terms. In particular, the cross term of disjoint signals (i.e. auto-terms do not overlap) has zero energy content. Figure 3.2 illustrates the energy content of the example above, and no energy exists for the cross term at t = 0. In this simple case the cross terms are separated spatially from the auto-terms. We can easily identify them and their integration in time has no contribution to total energy. However, for multicomponent or continuous frequency signals, cross-terms of their Wigner transform not only appear in the separate area, where the auto-terms are zero or negligible but also spread throughout the domain of the auto-terms, and corrupt them in the complete Wigner transform. Therefore, it is difficult to distinguish between auto-terms and cross terms if more complicated signals are used [START_REF] Hlawatsch | Linear and quadratic time-frequency signal representations[END_REF][START_REF] Ping | Cross-terms suppression in wigner-ville distribution based on image processing[END_REF].

W [F 1 ] in fig. 3.1b.
What calls for special attention is that a sum of complex-valued signal (3.11) was given above as an example. Solutions of wave propagation problems are always real-valued. It is known that any real-valued signal can be expressed as a sum of two complex-valued signals with positive and negative frequencies. For example, we have: cos(ωt) = e iωt +e -iωt 2 . Thus the cross terms in the Wigner transform of a sum of two real-valued signals result from not only interactions between intrinsic positive frequencies but also between positive and negative frequencies.

We study now the sum of two real-valued signals G = G 1 + G 2 given by:

G 1 (t) = F 1 (t) + F * 1 (t) 2 = A 1 h(t -t 1 ) cos(ω 1 t), G 2 (t) = F 2 (t) + F * 2 (t) 2 = A 2 h(t -t 2 ) cos(ω 2 t) (3.16)
with F 1 and F 2 already defined in eq. (3.11). So their complex conjugates F * 1 and F * 2 are respectively localized around (t 1 ; -ω 1 ) and (t 2 ; -ω 2 ) in the (t; ω)-plane. Thus G(t) is composed of four frequency components, and there are four auto-terms at (t 1 ; ±ω 1 ), (t 2 ; ±ω 2 ) and six cross terms at (t 1 ; 0), (t 2 ; 0), (t m ; ±ω m ), (t m ; ±ω d ) for its Wigner transform.

Considering the example of a sum of the two real-valued gaussian signals defined by eq. (3.16) with parameters given by eq. (3.15), its Wigner transform is computed and illustrated in fig. 3.3. We get four auto-terms at (t; ω) = (-2; ±20) and (2; ±40) and six cross terms at (t; ω) = (-2; 0), (2; 0), (0; ±30) and (0; ±10), which appear between each two auto-terms (fig. 3.3). Every cross term is resulting from the interference of each two auto-terms (marked with solid arrowed lines in fig. 3.3) and its locations, fluctuations and peak value are determined by eq. (3.13). Until now, we presented the existence and properties of cross terms in Wigner transform. Since the distribution of energy density in terms of frequency may be influenced by these terms, we will try and analyze some methods to remove them. However, the total energy property should not be influenced by that removal since in our work, the Wigner transform of wave fields is studied and the energy transport equation is used for evaluating errors.

For instance, cross terms between positive frequencies and negative frequencies for 1D signal can be suppressed by Hilbert transform. In fact, Hilbert transform of a signal allows eliminating its negative frequency components and we obtain what is called the analytical signal [START_REF] Hahn | Hilbert transforms in signal processing[END_REF][START_REF] Turner | Hilbert transforms, analytic functions, and analytic signals[END_REF]. Thus its Wigner transform contains only cross terms between non negative frequencies. The 1D form of analytical signal using the Hilbert transform has been used extensively in signal processing since their introduction by Gabor [START_REF] Gabor | Theory of communication. part 1: The analysis of information[END_REF]. However, their extension to the 2D case and their application to images are limited, due to the non-uniqueness of the multidimensional Hilbert transform. Some usual definitions of 2D Hilbert transform can be found in [START_REF] Lorenzo-Ginori | An approach to the 2d hilbert transform for image processing applications[END_REF]. However, either negative frequencies are not suppressed totally, or the real part of the constructed analytical signal does not coincide with the original signal. Consequently, since we are concerned with 2D or higher dimension, this method is not practical. Some suppression methods have been already developed to reduce the cross terms between positive frequencies, such as low pass filtering or adding window functions [START_REF] Pikula | A new method for interference reduction in the smoothed pseudo wigner-ville distribution[END_REF][START_REF] Shin | Pseudo wigner-ville time-frequency distribution and its application to machinery condition monitoring[END_REF], image processing [START_REF] Arce | Elimination of interference terms of the discrete wigner distribution using nonlinear filtering[END_REF][START_REF] Ping | Cross-terms suppression in wigner-ville distribution based on image processing[END_REF], combining the STFT and the Wigner transform [START_REF] Gonçalves | Pseudo affine wigner distributions: Definition and kernel formulation[END_REF][START_REF] Kadambe | A comparison of the existence of'cross terms' in the wigner distribution and the squared magnitude of the wavelet transform and the short-time fourier transform[END_REF][START_REF] Stanković | Time-frequency signal analysis based on the windowed fractional fourier transform[END_REF], etc. Generally, in all cases, cross terms can never be totally removed without affecting the appearance and quality of the resulting time-frequency representation [START_REF] Zou | Nonexistence of cross-term free time-frequency distribution with concentration of wigner-ville distribution[END_REF]. In detail, a loss of frequency concentration or resolution will appear with these methods and the most important energy properties cannot hold any more. Also, such methods are either computationally expensive or have a very limited applicability. Consequently, we do not use the classical cross term removal methods in our work so as to conserve always the energy properties. However, we should always pay attention to the existence of cross terms and their influence on our numerical computations, especially when explaining phenomena in phase space.

Scaled spatio-temporal Wigner transform in the high frequency limit

As discussed in the introduction of this chapter, we consider high frequency limit where the typical distance of propagation L of waves is much larger that the typical wavelength λ in the system:

ε = λ L ≪ 1.
The spatio-temporal Wigner transform (definition 3.1) should be rescaled and defined as follows:

Definition 3.3 (Scaled spatio-temporal Wigner transform). The scaled spatio-temporal Wigner transform W ε [u ε , v ε ] for two vector fields u ε (t, x) and v ε (t, x) is defined as:

W ε [u ε , v ε ](t, x; ω, k) := 1 ε 2 W[u ε , v ε ] t, x; ω ε , k ε = 1 (2π) 1+d R×R d u ε t - ετ 2 , x - εy 2 ⊗ v * ε t + ετ 2 , x + εy 2 e i(τω+y•k) dτdy (3.17)
and the (auto-) spatio-temporal Wigner transform

W ε [u ε ] := W ε [u ε , u ε ].
With this definition, the spatio-temporal Wigner transform is properly scaled at ε -1 for ω and k so as to observe all the fast fluctuations of order ε -1 and to separate them from slow fluctuations of ε 0 . It satisfies always all the properties of the temporal Wigner transform introduced in the section 3.1.1 and section 3.1.2.

Using rules of pseudo-differential calculus presented in Appendix C, we derive the following properties of the scaled Wigner transform: Energy-based error estimation for HF wave in homogeneous media Lemma 3.4 (Time and space derivative properties of Wigner transform).

∂ t W ε [u ε ] = W ε [∂ t u ε , u ε ] + (W ε [∂ t u ε , u ε ]) * T (3.18a) k • ∇ x W ε [u ε ] = W ε k • ∇ x u ε , u ε + W ε k • ∇ x u ε , u ε * T (3.18b)
where ( k • ∇ x W) i j = k l ∂ x l W i j (with Einstein summation convention). It is in fact the directional derivative along k of W.

Proof.

Using homogeneous pseudo-differential operator's properties eq. (C.4) and eq. (C.3) with ϕ(εD t ) = ε∂ t and ϕ(εD

x ) = ε k • ∇ x , we have    W ε [ε∂ t u ε , u ε ] = iωW ε [u ε ] + ε∂ t 2 W ε [u ε ] W ε [u ε , ε∂ t u ε ] = -iωW ε [u ε ] + ε∂ t 2 W ε [u ε ]    W ε ε k • ∇ x u ε , u ε = ikW ε [u ε ] + ε k•∇ x 2 W ε [u ε ] W ε u ε , ε k • ∇ x u ε = -ikW ε [u ε ] + ε k•∇ x 2 W ε [u ε ] (3.19)
Then taking the sum of each two equations in eq. (3.19), we get:

∂ t W ε [u ε ] = W ε [∂ t u ε , u ε ] + W ε [u ε , ∂ t u ε ] (3.20a) k • ∇ x W ε [u ε ] = W ε [ k • ∇ x u ε , u ε ] + W ε [u ε , k • ∇ x u ε ] (3.20b)
Recalling the conjugation property (eq. (3.4)) of Wigner transform), thus we obtain eq. (3.18). Now we present the Wigner measure, i.e. the energy density described by the transport equation or the radiative transfer equation. Definition 3.5 (Wigner measure). The Wigner measure, denoted by W, is the weak limit of the Wigner transform W ε as ε → 0:

W[u ε ] := lim ε→0 W ε [u ε ] (3.21)
As we observed with the examples in the section 3.1.2, Wigner transform can be positive or negative. However, Wigner measure can be demonstrated to be always positive [START_REF] Lions | Sur les mesures de wigner[END_REF]. Thus Wigner measure can be interpreted as the energy density of waves in phase space. This is exactly the property that makes the Wigner measure a useful tool for the study of wave or quantum interference phenomena. Furthermore, the potential and kinetic energy of high frequency propagating waves can be estimated by its Wigner measure with the following lemma: Lemma 3.6 (Energy property of Wigner measure of wave fields).

lim ε→0 1 2 R d (C : ε ε ε(u ε )(t, x)) : ε ε ε * (u ε )(t, x)dx = 1 2 R×R d ×R d ρΓ Γ Γ(k) : W[u ε ](t, x; ω, k)dωdxdk (3.22a) lim ε→0 1 2 R d ρ∂ t u ε (t, x) • ∂ t u * ε (t, x)dx = 1 2 R×R d ×R d ρTr (W[∂ t u ε ](t, x; ω, k)) dωdkdx (3.22b)
Proof. For eq. (3.22a), by using the integration property of Wigner transform over (ω, k) (the second equation of eq. (3.6)), we have

lim ε→0 1 2 R d (C : ε ε ε(u ε )(t, x)) : ε ε ε * (u ε )(t, x)dx = lim ε→0 1 2 R ((ε∇ x ) * • (C : ε ε ε(u ε )(t, x))) • u * ε (t, x)dx = lim ε→0 1 2 R×R d ×R d Tr (W ε [(ε∇ x ) * • (C : ε ε ε(u ε )(t, x)), u ε ] (t, x; ω, k)) dωdkdx (3.23)
Then applying the property of pseudo-differential operators eq. (C.3) by taking ψ(εD

x )u ε = (ε∇ x ) * • (C : ε ε ε(u ε ), we can derive: lim ε→0 1 2 R d (C : ε ε ε(u ε )(t, x)) : ε ε ε * (u ε )(t, x)dx = lim ε→0 1 2 R×R d ×R d ρ ρ -1 C : (ik + ε∇ x 2 ) * ⊗ (ik + ε∇ x 2 ) : W ε [u ε ](t, x; ω, k)dωdkdx = 1 2 R×R d ×R d ρΓ Γ Γ(k) : W[u ε ](t, x; ω, k)dωdkdx (3.24)
The kinetic energy eq. (3.22b) can be obtained by using directly the energy property of Wigner transform (the second equation of eq. (3.7)).

Note that these results in eq. (3.22a) and eq. (3.22b) are valid also in heterogeneous media when the elastic parameters are dependant of x [START_REF] Baydoun | Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media[END_REF].

In high frequency limit, highly oscillating quantities such as u 2 ε , σ σ σ 2 ε do not satisfy a closed-form equation in time-space (t, x); however, we will see that the Wigner measure W[u ε ] as ε → 0 satisfies a closed-form equation that we derive in next section 3.2.

Discrete Wigner transform and analysis of parameters

Before the introduction of the transport equation, numerical computation of Wigner transform is an important aspect in numerical applications. The Wigner transform of numerical solutions of wave equation is quantified in our work and its analytical expression does not exist. Besides, considering Energy-based error estimation for HF wave in homogeneous media the computational cost, Wigner transform should not be computed in an infinite region. Thus we should discuss the numerical computation of Wigner transform.

In this section, the discrete Wigner transform (DWT) is defined and analyzed. DWT for a temporal signal is presented here for simplicity (DWT for a time-space dependent function is realized in the same way). Numerical examples are given for Wigner transform of a ricker signal in time and Wigner transform of wave fields in the time-space domain.

The discretization of the temporal Wigner transform (3.2) is defined in the following: Definition 3.7 (Discrete Wigner transform). Assuming that a discretized signal is given as

F[N] := F(t = N∆t), N ∈ [1, N s ],
where N s is the total number of samples. ∆t is the sampling interval.

W[F][N; L] := W[F] t = N∆t; ω = 2π L N s ∆t = 1 π ∑ M F[N -M]F * [N + M]e i4πM L Ns (3.25)
Compared with the classical discrete Fourier transform (DFT) [START_REF] Bracewell | The Fourier transform and its applications[END_REF]:

G[L] = ∑ M G[M]e i2πM L Ns (3.26)
the DWT can be interpreted as a DFT of the correlation function FF * with twice frequency. Now we discuss the choice of numerical parameters in DWT: the sampling interval ∆t and the sampling length or computation length ζ t .

According to the Shannon sampling theorem [START_REF] Jerri | The shannon sampling theorem-its various extensions and applications: A tutorial review[END_REF], the maximum resolvable frequency must be half the sampling frequency f s for DFT. Thus for DWT, f s or ∆t satisfies:

f s = 1 ∆t ≥ 4 f F (3.27)
where f F is the frequency of signal F(t). So ∆t takes a very small value for high frequency f F .

Frequency resolution is defined as the distance between two adjacent frequency data points in the DWT so it stands for the precision of the frequency representation of the signal: Smaller ∆ f means a better frequency resolution. However, since ∆t is very small, increasing ζ t means increasing a large total number of sampling N s , that may cause a high computational cost especially for high dimension input. We have to choose an appropriate ζ t .

∆ f = 1 2ζ t = 1 2N s ∆t , f ∈ - 1 4∆t , 1
We consider firstly identifying the value of ζ t for two adjacent frequencies of a signal. If we want to separate two peaks with frequencies ω 1 and ω 2 of DWT (ω 2 > ω 1 ), we need to have at least one node between these two points in ω axis, i.e.:

ω 2 -ω 1 2π ≥ 2∆ f = 1 ζ t , i.e. ζ t ≥ 2π ω 2 -ω 1 (3.29)
Now considering a continuous frequency band with a maximum frequency ω 0 , it is known that it has a cross term ω 0 2 between ω 0 and low frequency near to zero in Wigner transform. For separating components ω 0 and ω 0 2 , we need

ζ t ≥ 4π ω 0 (3.30)
In conclusion, discrete Wigner transform in time is realized within a chosen computation length ζ t that satisfies eq. (3.30). Now same analysis of parameters can be extended to the discrete Wigner transform of a time-space dependant signal, denoted by (ζ t , ζ x ). They should verify:

ζ t ≥ 4π ω 0 , ζ x ≥ 4π k 0 (3.31)
where k 0 is the frequency related to the space variable, i.e the wave number for wave propagation problems. Since spatio-temporal Wigner transform depends on four parameters (t, x; ω, k), we fix one time-space point (t 0 , x 0 ) in time-space and represent all the results in phase space (ω, k), computed for the chosen computation region (ζ t , ζ x ) in time-space. In addition, we know that Wigner transform Energy-based error estimation for HF wave in homogeneous media is real (for d = 1) and symmetric about the origin in phase space so we illustrate only positive ω in all the figures. Now the DWT of a ricker signal is given so as to observe the influence of different computation lengths ζ t . It can be seen as the Wigner transform of the product of the ricker signal and a rectangular window function with different window lengths ζ t . Figure 3.4 illustrates a ricker signal centred at t = 0 s with period T R = 4 ms. The maximum frequency f max = 500 Hz (frequency for maximum amplitude). According to eq. (3.30), ζ t ≥ 2/ f max = 2T max = 4 ms for the identification of f max . Note also that 2T max = T R so this choice is reasonable in the sense that this lower bound allows including the whole signal F(t). Note also that since this ricker signal has a spectrum of frequency over a continuous frequency range, cross terms of its Wigner transform spreads out over the same continuous frequency range. Thus the maximum frequency found by its Wigner transform does not coincide with its real value f max = 500 Hz. To investigate the influence of the parameter ζ t on numerical calculations of Wigner transform, we propose to identify the maximum frequency by the middle point determined by the two adjacent frequencies corresponding to the two highest amplitudes (of course, the maximum at f = 0 due to cross terms with negative frequencies is neglected). It is respectively 375 Hz, 564 Hz, 530 Hz, 515 Hz, obtained with fig. 3.5f and illustrated in fig. 3.6 with ζ t increasing. It shows that the identified maximal frequency of the numerically calculated Wigner transform converges to a value that would be slightly larger than the maximal frequency f max defined by the ricker signal. Now to study the influence of the size of sampling length on the discrete spatio-temporal Wigner transform, we consider the wave propagation in a 1D semi-infinite homogeneous bar (fig. 3.7). Its end is subjected to Neumann boundary condition by a ricker signal F(t) with a magnitude of 1 kN/m 2 (t ∈ [0, T R ]). The related parameters are given in the table 3.1. Here the analytical solution of wave displacements can be computed with eq. (2.48).

-ζ t 2 , ζ t 2 × -ζx 2 , ζx 2 W[u](t 0 , x 0 ; ω, k)dωdk = |u(t 0 , x 0 )| 2 (3.32)
is always verified numerically, independent of the choice of computation region. It can be observed that in phase space, the spatio-temporal Wigner transform gives the energies along the direction of propagation (ω = c 0 |k|) in all cases.

Wigner transform applied to the elastic wave equation 57

In 1D homogeneous media, wave displacement is constant along the characteristics x = ±c 0 t. It is known that the Fourier transform of a constant function is a delta dirac function. It means that Wigner transform of wave front should be analytically a dirac function in the direction perpendicular to ω = -c 0 k. However, the numerical computation of Wigner transform along x = c 0 t is realized in a finite region. It can be seen as a Fourier transform for a rectangular function with different widths, so a sinc function is obtained in the direction perpendicular to ω = -c 0 k in all cases (fig. 3.8d).

Comparing the figures fig. 3.8a to fig. 3.8c, it is shown once again that the parameters (ζ t , ζ x ) = (8T max , 8λ max ) give a higher resolution of Wigner transform. Figure 3.8d gives the value of Wigner transform in the line perpendicular to ω = -c 0 k across the maximum frequency k = -k max = -157m -1 (black solid line in figs. 3.8a to 3.8c). It can be observed that the main lobe is more and more narrow, i.e. a better concentration, with the increasing of the size of computation region. As the computational cost with this choice is acceptable for our applications, with about 0.08 s for the calculation for one time-space point, thus these parameters are used for all the calculations presented in the section 3.3 and in the chapter 4.

Transport energy-based residual errors of numerical solutions of wave equation

In this section, the transport equation with the scaled Wigner transform is derived in the weak coupling regime, based on the wave equation. It reveals the connection of the wave equation to the energy transport equation in homogeneous media [START_REF] Bal | Kinetics of scalar wave fields in random media[END_REF][START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF]. Then residual-type errors of finite element solutions of the wave equation are defined based on the transport equation in terms of Wigner transform of wave fields.

Transport equation in homogeneous media

The wave equation (1.1) in homogeneous media is considered here with assumption of f = 0. We apply the pseudo-differential calculus for this equation in the following.

Γ Γ Γ(ε∇ x ) denotes the pseudo-differential operator defined using the same mathematical symbol as the Christoffel operator Γ Γ Γ(iεk) (see Appendix C for the definition of the pseudo-differential operator and its symbol). Considering homogeneous and isotropic media, we have:

Γ Γ Γ( k) = c 2 s I d + (c 2 p -c 2 s ) k ⊗ k (3.33)
We rescale t → t/ε and x → x/ε in high frequency limit. Then the wave equation can be expressed as:

(ε∂ t ) 2 I d -Γ Γ Γ(ε∇ x ) • u ε (t, x) = 0 (3.34)
and the cross Wigner transform is applied to this equation and u ε :

W ε (ε∂ t ) 2 I d -Γ Γ Γ(ε∇ x ) • u ε , u ε = 0, W ε u ε , (ε∂ t ) 2 I d -Γ Γ Γ(ε∇ x ) • u ε = 0 (3.35)
Using the properties of pseudo-differential operators in eq. (C.3), eq. (C.4), and eq. (C.2) by taking ϕ(εD t ) = (ε∂ t ) 2 and ψ(εD x ) = Γ Γ Γ(ε∇ x ), we obtain:

(iω) 2 W ε [u ε ] + iεω∂ t W ε [u ε ] = Γ Γ Γ (ik) • W ε [u ε ] - iε 2 {Γ Γ Γ (ik) , W ε [u ε ]} + o(ε 2 ), (iω) 2 W ε [u ε ] -iεω∂ t W ε [u ε ] = W ε [u ε ] • Γ Γ Γ ((ik)) * T + iε 2 W ε [u ε ], (Γ Γ Γ (ik)) * T + o(ε 2 ) (3.36)
where the Poisson bracket is defined by:

{Γ Γ Γ, W} = ∇ k Γ Γ Γ • ∇ x W -∇ x Γ Γ Γ • ∇ k W with the product ∇ k Γ Γ Γ • ∇ x W := ∂ k j Γ Γ Γ • ∂ x j W (with Einstein summation convention).
The asymptotic expansion of W ε is defined by:

W ε := W 0 + εW 1 + o(ε 2 ) (3.37)

Transport energy-based residual errors of numerical solutions of wave equation 59

Note that the 0-order term W 0 is equivalent to the Wigner measure W, as W 0 = lim ε→0 W ε . Introducing eq. (3.37) into eq. (3.36), the leading order o(ε 0 ) gives:

ω 2 I d -Γ Γ Γ(k) W 0 = 0 (3.38)
Using the eigendecomposition for Γ Γ Γ as we did for the Christoffel equation (1.8), it can be found that the non-trivial solutions of eq. (3.38) can be written in the form [START_REF] Bal | Kinetics of scalar wave fields in random media[END_REF]:

W 0 = W p k ⊗ k + W s = W p k ⊗ k + ∑ α,β =1,2 W s,αβ k⊥ α ⊗ k⊥ β (3.39)
and the eigenvalues of the Christoffel tensor:

ω 2 p = c 2 p |k| 2 , or ω 2 s = c 2 s |k| 2 (3.40)
In eq. (3.39), W 0 is decomposed into P-mode (associated to the eigenvalue ω p and the eigenvector k) and S-mode (associated to the multiple eigenvalue ω s and the eigenvectors k⊥ 1 , k⊥ 2 ). W p is a scalar function and W s is a 2 × 2 matrix (two modes for S-wave) in 3D media. Thus the multiply-scattered wave energy in an elastic medium may be characterized by four parameters, one for the longitudinal wave and three for the transverse wave (according to the symmetry property of Wigner transform, only three coefficients of W s are independent). Now by subtracting the two equations in eq. (3.36) and inserting eq. (3.37), the order o(ε) of the obtained equation gives:

2ω∂ t W 0 = 1 2 Γ Γ Γ (k) , W 0 -W 0 , Γ Γ Γ (k) (3.41)
We multiply the above equation on the left and the right by k for mode P and by k⊥ for mode S, and recall that the identity matrix

I d = k ⊗ k + I ⊥ = k ⊗ k + k⊥ 1 ⊗ k⊥ 1 + k⊥ 2 ⊗ k⊥ 2 .
Its right-hand side reads:

                           k • Γ Γ Γ (k) , W 0 + W 0 , Γ Γ Γ (k) • k = k • Γ Γ Γ (k) • I d , W 0 • k -k • W 0 , I d • Γ Γ Γ (k) • k = ω 2 p , W p -W p , ω 2 p = 2∇ k ω 2 p • ∇ x W p I ⊥ • Γ Γ Γ (k) , W 0 + W 0 , Γ Γ Γ (k) • I ⊥ = I ⊥ • Γ Γ Γ (k) • I d , W 0 • I ⊥ -I ⊥ • W 0 , I d • Γ Γ Γ (k) • I ⊥ = ω 2 s , W s -W s , ω 2 s = 2∇ k ω 2 s • ∇ x W s (3.
42) Then introducing eq. (3.42) into eq. (3.41) and using the definition of ω p and ω s , the following equations for mode P and S can be found:

   ω∂ t W p [u ε ] -c 2 p k • ∇ x W p [u ε ] = 0 ω∂ t W s [u ε ] -c 2 s k • ∇ x W s [u ε ] = 0 (3.43)

Energy-based error estimation for HF wave in homogeneous media

We recall that k • ∇ x is the operator of the directional derivative along k. According to eq. (3.40), for example in the case of P-waves, we have ω p = ±c p |k|. Considering that ω is physically always positive and using the symmetry about zero of Wigner transform in phase space, we define ω > 0 and ±k that represent the two directions "±". Thus W p [u ε ] can be projected into these two directions by:

W p [u ε ](t, x; ω, k) = W + p [u ε ](t, x; -k)δ (ω -c p (-k • k)) + W - p [u ε ](t, x; k)δ (ω -c p (k • k)) (3.44)
This shows that the forward wave intensity W + p propagates in the direction -k and the backward wave intensity W - p propagates in the direction k in phase space. Note that this definition is consistent with the numerical example in the section 3.1.4. In fig. 3.8, it is noticed that Wigner transform of forward waves (t, x = c p t > 0) in 1D media takes value with slope of -c p in phase space (for k < 0, ω = -c p k > 0). In fact, it is a result of exponential terms e j(τω+yk) in the definition of Wigner transform, which differs with the common definition of harmonic forward-going plane waves e j(τω-yk) . However, it does not affect the analysis of results. We keep this traditional definition of Wigner transform and use this convention of sign for (ω, k) in eq. (3.44) in the following. All analysis and definitions above are applied similarly for S-wave.

Introducing eq. (3.44) into eq. (3.43) (note that k in eq. (3.43) should be changed to -k for W + p and W + s according to our choice of sign), it can be decomposed into:

   ∂ t W ± p [u ε ] ± c p k • ∇ x W ± p [u ε ] = 0 ∂ t W ± s [u ε ] ± c s k • ∇ x W ± s [u ε ] = 0 (3.45)
which are the transport equations for P and S-wave in homogeneous media, subjected to some initial conditions and boundary conditions.

The general solutions of eq. (3.45) for harmonic waves can be expressed as:

   W ± p (t, x; ∓k) = A p ( k • x ∓ c p t) W ± s (t, x; ∓k) = A s ( k • x ∓ c s t) (3.46)
In the time-space domain (t, x), general solutions of P-wave displacement in homogeneous media

∂ 2 t u -c 2 p ∇ 2
x u = 0 can be expressed as:

u p (t, x) = F( k • x -c p t) + G( k • x + c p t) (3.47)
where F and G are arbitrary functions of their argument. It has two components that represent respectively the propagation of wave fronts in the forward and backward directions. It can be remarked that the energy density W p (eq. (3.46)) has the same arguments in the time-space domain.

Considering the Wigner transform of u p (t, x), according to eq. (3.10), it has two auto-terms related to energy density of two directions W ± p in eq. (3.46) and one cross term that vanishes in the weak sense when ε → 0. The same can be said of the S-wave. In the phase space (ω, k), Wigner measure presents the relation of dispersion along the characteristics in eq. (3.44), that is verified numerically by the example in fig. 3.8 of section 3.1.4.

Theoretical residual errors based on transport equation

We assume that u h∆t is a finite element solution of the wave equation (1.1) in the high frequency limit. The numerical errors related to u h∆t are evaluated based on transport equations eq. (3.45) in terms of energy quantities of wave fields, i.e. the projections of Wigner transform of u h∆t . In the following, its projection on mode P is denoted by W

p [u h∆t ] = k • W ε [u h∆t ] • k and on mode S by W s [u h∆t ] = I ⊥ • W ε [u h∆t ] • I ⊥ .
Definition 3.8. Residual errors based on transport equations (3.45) are defined as:

   R ± p [u h∆t ](t, x; ω, k) := ∂ t W ± p [u h∆t ](t, x; ω, k) ± c p kp • ∇ x W ± p [u h∆t ](t, x; ω, k) R ± s [u h∆t ](t, x; ω, k) := ∂ t W ± s [u h∆t ](t, x; ω, k) ± c s ks • ∇ x W ± s [u h∆t ](t, x; ω, k) (3.48)
Furthermore, it can be demonstrated that eq. (3.48) can be transformed into the following theorem.

Theorem 3.9. The local residual error maps in (t, x; ω, k) domain can be obtained by:

   R ± p [u h∆t ] = W ± p [v h∆t , u h∆t ] + W ± p [v h∆t , u h∆t ] * T ± c p W ± p [ kp • ∇ x u h∆t , u h∆t ] + W ± p [ kp • ∇ x u h∆t , u h∆t ] * T R ± s [u h∆t ] = W ± s [v h∆t , u h∆t ] + (W ± s [v h∆t , u h∆t ]) * T ± c s W ± s [ ks • ∇ x u h∆t , u h∆t ] + W ± s [ ks • ∇ x u h∆t , u h∆t ] * T (3.49)
Proof. Using lemma 3.4.

Errors defined with eq. (3.49) is more interesting than eq. (3.48) because it emphasizes the importance of the correlations of the displacement with the velocity and the directional derivative along the wave vector of the displacement. Numerically, it avoids the computation of the time derivatives by finite difference, when numerical velocities are already obtained at each time step.

Considering that the local errors R ± [u h∆t ](t, x; ω, k) depend on four parameters and using the energy properties of Wigner transform (eq. (3.7)), we propose two new error quantities in an energy norm:

R[u h∆t ](ω, k) := R×R d ∑ α tr R + α [u h∆t ](t, x; ω, k) + R - α [u h∆t ](t, x; ω, k) dtdx (3.50) and R[u h∆t ](t, x) := R×R d ∑ α tr R + α [u h∆t ](t, x; ω, k) + R - α [u h∆t ](t, x; ω, k) dωdk (3.51)
where α = P, S stands for different modes. Taking into account that the existence of cross terms disturbs the distribution of energy in terms of frequency, the second definition R(t, x) is used and computed in this work.

Energy-based error estimation for HF wave in homogeneous media

Numerical application for wave propagation in a 1D homogeneous medium

In this section, we analyze numerically the defined residual errors eq. (3.54) and eq. (3.55) in a 1D case.

We firstly rewrite the transport equation and the residual errors in 1D media for clarity and simplification of some notations. In this case, only P-wave propagates with constant velocity c 0 and the subscript p can be neglected. The transport equation eq. (3.45) reads:

∂ t W ± [u ε ] ± c 0 ∂ x W ± [u ε ] = 0 (3.52a)
Recall that W ± with the convention of sign defined in eq. (3.44) can be written as:

W + [u ε ] := W[u ε ](t, x; ω = -c 0 k > 0, k < 0), W -[u ε ] := W[u ε ](t, x; ω = c 0 k > 0, k > 0) (3.53)
The residual error maps in 1D are obtained with eq. (3.49) and eq. (3.51), namely:

R ± [u h∆t ](t, x; ω, k) := 2Re W ± ε [v h∆t , u h∆t ](t, x; ω, k) ± 2c 0 Re W ± ε [ε h∆t , u h∆t ](t, x; ω, k) (3.54) and R[u h∆t ](t, x) := R×R R + [u h∆t ](t, x; ω, k) + R -[u h∆t ](t, x; ω, k) dωdk (3.55)
where v h∆t and ε h∆t stand for respectively the numerical velocity and strain fields. Now we consider here the same 1D wave propagation model as the one defined in the section 3.1.4 (fig. 3.7). All the physical parameters are defined already in the table 3.1 except that here we consider a bar with a finite length L = 1 m so wave front propagates a round trip in a time period T f = 0.4 ms. The left end is subjected to the Neumann conditions with the ricker signal and the right end to the free boundary condition.

The exact solutions of the wave displacement in the time-space domain are obtained analytically and illustrated in fig. 3.9. Note that the Wigner transform of the exact solutions is calculated numerically.

In order to analyze the evolution of residual errors with refinement, three numerical finite element solutions are computed with different mesh sizes h and time steps ∆t shown in the table 3.2. Note that they all satisfy the sampling condition in eq. (3.27): 1 ∆t ≥ 4 f c , with f c the cutoff frequency of the ricker signal, which can be chosen equal to 3 f max . Generally, the maximum of the estimated errors is located around k max (we ignore the maximum at k = 0 where all the cross terms between positive/negative frequencies overlap). In the case λ ( f c )/h = 4 for point B (the top figure in fig. 3.11b), it is observed that the frequency band is much shorter than other cases and the maximum of errors is approximatively located at k max /2. In fact, it can be remarked in fig. 3.12 (more evident in fig. 3.12b) significant errors related to the numerical diffusion of numerical solutions with the time discontinuous Galerkin method used here. It includes the dissipation (associated to the amplitude error) and the dispersion (associated to the phase error). These effects, depending on mesh size and the size of the time step, increase with time. Consequently, for the coarse refinement (λ ( f c )/h = 4 here), the quality of estimated errors cannot be ensured. Errors can be underestimated (fig. 3.13b). More studies on diffusion properties of DG or other numerical methods for the wave equation can be found in [START_REF] Ainsworth | Dispersive and dissipative properties of discontinuous galerkin finite element methods for the second-order wave equation[END_REF][START_REF] Hu | An analysis of the discontinuous galerkin method for wave propagation problems[END_REF][START_REF] Steffens | A simple strategy to assess the error in the numerical wave number of the finite element solution of the helmholtz equation[END_REF]. Some noticeable fluctuations are also found along ω = ∓c 0 k (most remarkable for λ ( f c )/h = 4). At the bottom of fig. 3.13, we give respectively the interpolation of the previous three plots on the line perpendicular to ω = ∓c 0 k across the k = ∓k max for point A and B. For point A, in the case of λ ( f c )/h = 4, it is observed that these fluctuations are small compared to the relatively large errors in ω = c 0 k. In the case of λ ( f c )/h = 20, fluctuations and errors become more comparable. These fluctuations are resulting from "spectral leakage". In simple terms, DFT (so does DWT) for the sampled signal is repeated periodically; if some glitches or discontinuities exist in the boundaries of the chosen computation region, a leakage of frequency will arise. Leakage causes the signal levels to be reduced and redistributed over a broad frequency range. For exact solutions of wave propagation in homogeneous media, if we compute its Wigner transform along the direction of propagation x = ±c 0 t, it can be seen as a constant continuous function. No leakage happens. For numerical solutions, discontinuities always exist and we cannot find a method that allows removing totally these fluctuations. Increasing the size of computation region is one of methods to reduce them. That is also an explanation for choosing a large (ζ t , ζ x ) = (8T max , 8λ max ). A detailed discussion about this phenomenon can be found in appendix B.

Numerical application for wave propagation in a 1D homogeneous medium

63 (t A ,x A ) ζ x ζ t (t B ,x B ) ζ x ζ t (a)

Residual errors in time-space for different refinements of finite element solutions

Considering the numerical fluctuations found in last section, R(t, x; ω, k) is filtered with a window function and then we compute errors in the time-space domain R(t, x) defined in eq. (3.55), i.e.

R(t, x)

:= R×R R ± (t, x; ω, k) h R (ω, k)dωdk (3.56)
where h R (ω, k) is a window function defined in phase space, which allows the filtering of fluctuations according to our choice. Considering here that the integral area of errors should be only around ω = ±c 0 k so we choose naturally a rectangular window defined with the rotational coordinate (it can be also found in Appendix F for a more detailed definition of this new coordinate):

hR (k ⊥ , k ∥ ) = h R (ω, k) = rect(k ⊥ )rect(k ∥ ) with    k ⊥ = ω c 0 -k √ 2 k ∥ = ω c 0 +k √ 2 (3.57)
For forward waves (point A, fig. 3.14), the length of hR along k ⊥ is from zero to the cutoff frequency:

k ⊥ ∈ [-3 √ 2k max , 0]. Its width in k ∥ direction is carefully chosen as ζ ∥ = 2 √ 2∆k
(recall that ∆k is frequency resolution). This width covers at least the main lobe of Wigner transform in this direction. It allows to contain the residual errors along the characteristic line but not too much fluctuations. Indeed, as we discussed in the section 3.1.3, the DWT can be seen as DFT of a rectangular function with period of (c

0 ζ t ) 2 + ζ 2 x = √ 2ζ
x in the (c 0 t, x)-plane along x = c 0 t (t axis is rescaled to c 0 t) in homogeneous media, so the width of main lobe is equal to 2 Now residual errors are computed for all points (t, x) so we get a time-space maps of errors. Note that the boundaries are neglected since in the present work the chosen calculation domain of Wigner transform and the validity of the transport equation are only considered in the bulk of the medium. In fig. 3.15, the convergence of errors with decreasing mesh size is illustrated. It can be also found in fig. 3.15d that the integral of errors over space is increasing with time (more evident for the coarsest mesh). To understand more clearly, we compare values of all R(t, x) for a certain position x 0 or a certain moment t 0 in one plot. In fig. 3.16, it is observed that the residual errors defined in eq. (3.55) are decreasing with smaller mesh size as we expected. Two errors for exact solutions are given (discretized both with respect to the finest refinement λ /h = 20): one is computed in coordinates (t, x), another is along the direction of propagation. As we mentioned in last section, Wigner transform of the latter has no leakage. Thus it is noticed in fig. 3.16 that it has much smaller value than the former case and other numerical solutions. It means that the window function h R allows reducing the fluctuations but it cannot remove totally the effect of frequency leakage. We can conclude that errors resulting from leakage of frequency are very important in numerical computations of Wigner transform but it does not influence the evaluation of errors with refinement here. Finally, for only exact solutions without leakage in FFT, sources of residual errors in numerical computations can be considered mainly as the discretization in the time-space domain of analytical Energy-based error estimation for HF wave in homogeneous media equations. Indeed, errors of exact solutions with different sizes of discretization can be also computed and compared in fig. 3.17. We find that they are much smaller than errors in fig. 3.16 so they can be neglected.

Energy-based error estimation for HF wave in homogeneous media

(-3k max ,0) k || k ⊥ ζ | |

Conclusion

In this chapter, we presented the transport equation which describes the energy transport for HF wave propagation in homogeneous media, derived by the scaled spatio-temporal Wigner transform. The properties of cross terms in Wigner transform were analyzed and discussed in detail. For numerical computations of Wigner transform, an important study of parameters, i.e. the size of computation region in time and in space for discrete Wigner transform was determined in function of signal frequency considering the frequency resolution. Then the energy-based residual errors of the wave equation with the transport equation were proposed analytically. A numerical application of these results is realized in a 1D homogeneous medium. The errors were computed for approximate solutions with different mesh sizes and we compare all these results in the time-space domain. We observed that with increasing number of degrees of freedom, the defined residual errors are decreasing as we expected, which allows validating our tools of error estimator.

Chapter 4

Energy-based a posteriori error estimation for high frequency wave equation in heterogeneous media This chapter deals with high frequency wave propagation problems in heterogeneous media. Compared to the transport equation in chapter 3, a more general and complicated radiative transfer equation is studied. It is derived from the wave equation in the same way as we did in the section 3.2.1. It describes not only the transport of energy, but also the scattering of energy by heterogeneities, with the effects of coupling between different polarization modes. After a long propagation time, we find the mean amplitude of waves decreases with the distance travelled, since the coherent or mean wave energy is converted to incoherent fluctuations. Considering the difficulties and high computational costs of analyzing errors for highly oscillating wave fields in these media, we propose here to quantify Energy-based error estimation for HF wave in heterogeneous media errors of finite element solutions of the wave equation in terms of energy quantities governed by this new equation.

In the section 4.1, we start with the introduction of random media, which is used to model the heterogeneities, and the radiative transfer equation in these media. Then the analytical solutions of the radiative transfer equation in 1D media are obtained. The related propagation phenomena characterized by the scattering mean free length, are discussed. A numerical example of wave propagation in a 1D medium is presented and the localization phenomenon is discussed.

In the section 4.2, the energy-based residual errors of wave fields using the radiative transfer equation are defined analytically and computed numerically for FE wave solutions obtained in 1D heterogeneous media.

In the section 4.3, a filtering of the obtained residual errors with window functions is proposed, which allows smoothing Wigner transform and reducing the influence of cross terms. Its numerical application is presented and analyzed.

Finally, in the section 4.4, the energy-based residual errors of numerical wave solutions obtained with two different refinements are compared.

Radiative transfer equation in heterogeneous media in the weak coupling limit

In this section, we present the radiative transfer equation in heterogeneous media and the related phenomena. Firstly, the heterogeneities of media are described by a random model of mechanical properties. Then the scaled spatio-temporal Wigner transform defined in the chapter 3 is used for the derivation of the radiative transfer equation in heterogenous media. Some important characteristic lengths and propagation phenomena are discussed in the end.

Statistical description of weak randomly-fluctuating media

As we introduced in the section 1.2.4, we consider the high frequency wave propagation in the weak coupling limit. The physical properties of weakly heterogeneous media are assumed in the form:

C(x) = C 0 (x) + √ εC 1 x ε (4.1)
Here ρ(x) is assumed slowly varying so only the elastic tensor is assumed to vary rapidly.

The slowly varying part C 0 (x) satisfies:

C 0 (x) = E [C(x)] (4.2)
where E[•] denotes mathematical expectation or ensemble average.

The fast fluctuating part C 1 (y), with y = x ε , is considered as a statistically homogeneous mean zero random field with unit-variance in the y variable. The standard deviation of C is given by the mean square of the fractional fluctuations √ ε before C 1 . This size is the unique scaling which allows significantly modifying the energy spreading in the transport regime at a long propagation distance. Indeed, if the random fluctuations are too weak they will not affect significantly energy transport (although strong interactions between waves and media exist in this case with l c ≈ λ ), and too large fluctuations will lead to the localization of wave energy in media where the radiative transfer theory cannot be applied [START_REF] Sheng | Introduction to wave scattering, localization and mesoscopic phenomena[END_REF]. Note that the period of the fluctuations of the properties is of order ε, as the characteristic length over which the material fluctuates. This model allows an asymptotic analysis of wave propagation problems.

A random medium can be seen as an ensemble of heterogeneous media, called realizations. Each realization differs from another in the detailed spatial structure of fluctuations, but they have some common statistical properties, such as mean value, variance or standard deviation, and covariance function.

Definition 4.1. For any statistically homogeneous random field U, its covariance function P(r) and power spectral density P(k) (Fourier transform of covariance function) are respectively defined by:

P(r) = P(y -z) = E [U(y) ⊗ U(z)] , P(k) = R d P(z)e -iz•k dz (4.3)

Energy-based error estimation for HF wave in heterogeneous media

With the definition of P(k), we have:

(2π) d δ (p + q) P(p) = E U(p) ⊗ U(q) (4.4)
The covariance function is a statistical measure of the spatial correlation and the magnitude of the fluctuations in the medium. When a medium is isotropic, the covariance function depends only on the lag distance r = |y -z|. Three covariance functions are given here as examples: the sinc 2 , the exponential and the gaussian functions. These three functions and their corresponding power spectral densities are given in the table 4.1 and illustrated in fig. 4.2. Here "sinc" stands for the unnormalized cardinal sine function sin(x) x and "tri" stands for the scaled triangular function defined by tri(k/a) := 1 -|k| /a |k| ≤ a 0 otherwise . type covariance function power spectral density It can be seen that the covariance function P(r) depends on two parameters: the variance σ 2 and the correlation length l c . Here l c is defined so that it satisfies:

sinc 2 P(r) = σ 2 sinc 2 (πr/l c ) P(k) = σ 2 l c tri(kl c /2) exp P(r) = σ 2 e -2r/l c P(k) = 1 π σ 2 l c 1 (1+k 2 l 2 c /4) 2 gauss P(r) = σ 2 e -πr 2 /l 2 c P(k) = 1 π σ 2 l c e -k 2
l c = 2 ∞ 0 P(r) σ 2 dr (4.5)
Different definitions for l c can be found in the literature [START_REF] Codona | Electromagnetic wave propagation through random media[END_REF][START_REF] Shinozuka | Response variability of stochastic finite element systems[END_REF]. Physically, the correlation length l c refers to an order of magnitude of the distance beyond which the values of the random field are almost uncorrelated.

Briefly, a random medium realization is obtained by the filtering of a white noise. This white noise is filtered with the square root of the power spectral density of the required covariance function in the wave number domain. The inverse Fourier transform provides the fluctuating field which is superposed on the slowly fluctuating background. Figure 4.2 gives one realization of numerical 2D random fields with spectra P(k x , k y ) = σ 2 l 2 c tri(k x l c /2)tri(k y l c /2) by spectral approaches [START_REF] Shinozuka | Simulation of stochastic processes by spectral representation[END_REF][START_REF] Shinozuka | Simulation of multi-dimensional gaussian stochastic fields by spectral representation[END_REF]. The influence of the correlation length on media can be observed.

Radiative transfer equation

In this section, the radiative transfer equation is derived, which describes the high frequency wave propagation in random media in the weak coupling limit. It provides a description of the evolution of wave energy in phase space in terms of the Wigner measure.

We present hereafter only the main steps taken in the derivation of the radiative transfer equation. More details can be found in [START_REF] Bal | Kinetics of scalar wave fields in random media[END_REF][START_REF] Baydoun | Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media[END_REF][START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF]. The main idea is the same as the derivation of the transport equation presented in the section 3.2.1. However, a special attention should be paid for the new order √ Energy-based error estimation for HF wave in heterogeneous media

Γ Γ Γ • u ε = ρ -1 (x) (C 0 (x) : (k ⊗ s u ε )) • k Γ Γ Γ 0 (x;k) + √ ερ -1 (x) (C 1 (x/ε) : (k ⊗ s u ε )) • k Γ Γ Γ 1 (x,x/ε;k) + ερ -1 (x) (∇ x C 0 (x) : (k ⊗ s u ε )) • k Γ Γ Γ 2 (x;k) +o(ε 3 2 ) (4.7)
The following equations hold when the spatio-temporal Wigner transform is applied to the wave equation and the wave field:

W ε [ (ε∂ t ) 2 I d -Γ Γ Γ(x; ε∇ x ) • u ε , u ε ] = 0, W ε [u ε , (ε∂ t ) 2 I d -Γ Γ Γ(x; ε∇ x ) • u ε ] = 0 (4.8)
A two-scale asymptotic expansion of the Wigner transform for these two equations is defined as:

W ε (t, x, x/ε; ω, k) = W 0 (t, x; ω, k) + √ εW 1 (t, x, x/ε; ω, k) + εW 2 (t, x, x/ε; ω, k) + o(ε 3
2 ) (4.9)

Introducing eq. (4.7) and eq. ( 4.8) into eq. ( 4.6), all the terms inside the Wigner transform and the Christoffel tensor are developed and expanded. Parameter identification is then applied respectively for the first three orders o(ε 0 ), o( √ ε) and o(ε).

The order o(ε 0 ) gives:

ω 2 I d -Γ Γ Γ 0 (x; k) • W 0 = 0 (4.10)
In the case of Γ Γ Γ 0 (x) is independent of x, the same equation as that in homogeneous media (eq. (3.38)) is found. In the general case, we also use the same eigendecomposition of W 0 in the P-mode and the S-mode as eq. (3.39), and the dispersion relation eq. (3.40) with c p = c p (x) and c s = c s (x).

The order o( √ ε) allows deriving the Fourier transform of W 1 with respect to y and the order o(ε) gives an equation in terms of W 0 and W 1 , while W 2 disappears under an assumption of orthogonality. Note that an ensemble average is applied when we deal with the order o(ε). The projection of W 0 on a mode α, denoted by W α (α = p, s), obtained in the radiative transfer equation is in fact E[W α ]. However, we denote it always by W α in the rest of document. Besides, a crucial mixing assumption is applied during the derivation: the average of the leading term W 0 is assumed to depend only on slow space variable x [START_REF] Bal | Kinetics of scalar wave fields in random media[END_REF][START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF]. It allows us to separate the product of the average of W 0 and the average of fast fluctuating variables since they vary on different scales. Finally combing equations in these two orders and using all these assumptions presented above, the 3D radiative transfer equation for the P-mode and the S-mode is obtained and given in the following:

∂ t W ± p (t, x; k) ± ω p (x; k), W ± p (t, x; k) = ω p (k)=ω p (k ′ ) s pp (x; k, k ′ )W ± p (t, x; k ′ )dk ′ -S pp (x; k)W ± p (t, x; k) + ω p (k)=ω s (k ′ ) s ps (x; k, k ′ ) k⊥ ′ ⊗ k⊥ ′ : W ± s (t, x; k ′ ) dk ′ -S ps (x; k)W ± p (t, x; k) (4.11a) ∂ t W ± s (t, x; k) ± ω s (x; k), W ± s (t, x; k) + N(x; k), W ± s (t, x; k) = ω s (k)=ω s (k ′ ) s ss (x; k, k ′ )G(k, k ′ ) • W ± s (t, x; k ′ ) • G(k, k ′ ) T dk ′ -S ss (x; k)W ± s (t, x; k) + ω s (k)=ω p (k ′ ) s sp (x; k, k ′ ) k′ ⊥ ⊗ k′ ⊥ W ± p (t, x; k ′ )dk ′ -S sp (x; k)W ± s (t, x; k) (4.11b) with G(k, k ′ ) = ( k • k′ ) k⊥ i ⊗ k⊥ i,⊥ ′ + k′ ⊥ ⊗ k⊥ ′ , i = 1, 2 (
with Einstein summation convention) (4.12)

S αβ (x; k) = ω α (k)=ω β (k ′ ) s αβ (x; k, k ′ )dk ′ , α, β = p, s (4.13)
Recall that ( k, k⊥ 1 , k⊥ 2 ) forms an orthonormal triplet. In the eqs. (4.11) and (4.12), the subscripts ⊥ and ⊥ ′ denote respectively the projection into the plane perpendicular to k and to k′ . It means:

k⊥ ′ = k -( k • k′ ) k′ , k′ ⊥ = k′ -( k′ • k) k, k⊥ i,⊥ ′ = k⊥ i -( k⊥ i • k′ ) k′ (4.14)
[A, B] := AB -BA stands for the Lie bracket. The matrix N(x; k) in eq. (4.11b) leads to the coupling of elements of the matrix W s . It is related to the slow variations of the background C 0 (in proportion to ∂ c s ∂ x j ) and vanishes in homogeneous media. s αβ (x; k, k ′ ) is the differential scattering cross-section for the mode αβ scattering, namely the rate of conversion of energy at x from the mode α with wave number k ′ into the mode β with wave number k ′ . It is associated to the fast fluctuations of media. S αβ (x; k) is the total scattering cross-section. For simplicity, in order to avoid the multiple S-mode in 3D media, the formula of differential scattering cross-sections in 2D isotropic media are given here. Results for 3D or anisotropic media can be found in [START_REF] Baydoun | Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media[END_REF]. The elasticity tensor C(x) is determined by two Lamé's coefficients λ and µ in the same form with eq. (4.1):

λ (x) = λ 0 (x) + √ ελ 1 x ε , µ(x) = µ 0 (x) + √ ε µ 1 x ε (4.15)
The differential scattering cross-sections are given by :

s pp (x; k, k ′ ) = c p (x) |k| 2 4 (c 2 p (x) -2c 2 s (x)) 2 c 4 p (x) P λ λ (k-k ′ )+ 4c 2 s (x)(c 2 p (x) -2c 2 s (x)) c 4 p (x) ( k• k′ ) 2 P λ µ (k-k ′ ) + 4c 2 s (x) c 4 p (x) ( k • k′ ) 4 P µ µ (k -k ′ ) (4.16) s ps (x; k, k ′ ) = c s (x) |k| 2 4 4( k • k′ ) 2 P µ µ (k -k ′ ) (4.17) s ss (x; k, k ′ ) = c s (x) |k| 2 4 P µ µ (k -k ′ ) (4.18) s sp (x; k, k ′ ) = c s (x) |k ′ | 2 4 4( k • k′ ) 2 P µ µ (k -k ′ ) (4.19)
where P λ λ , P µ µ , P λ µ are the power spectral densities of the fast fluctuations of corresponding subscript. For example, P λ µ is defined by (variance σ 2 = ε):

(2π) 2 δ (p + q) P λ µ (p) = εE λ 1 (p) µ 1 (q) (4.20)
Equation (4.16) implies that the fast variations C 1 contribute to the right-hand side of the radiative transfer equation only through the power spectral densities of its elasticity parameters.

Compared with the transport equation in homogeneous media (3.45), in the radiative transfer equations eq. ( 4.11) we have two additional terms on the left-hand side (transport of waves), {ω α , W ± α } (α = p, s) and [N, W ± s ], resulting from slow variations C 0 , and two additional terms on the right-hand side (scattering of waves), resulting from fast variations C 1 . Now the particular case of 1D heterogenous media with only fast fluctuations, is discussed. We assume that ρ and E 0 are constant, the heterogeneity in eq. ( 4.1) becomes:

E(x) = E 0 + √ εE 1 x ε (4.21)
The only mode in 1D is a P-mode with velocity c 0 = E 0 ρ , the radiative transfer equation eq. (4.11) becomes (the subscript α for mode P is removed hereafter for simplicity):

∂ t W ± (t, x; k) ± c 0 ∂ x W ± (t, x; k) = |k|=|k ′ | s(k, k ′ )W ± (t, x; k ′ )dk ′ -S(k)W ± (t, x; k) (4.22)
The differential scattering cross-section s(k, k ′ ) (4.16) and the total scattering cross-section S(k) (4.13) are:

s(k, k ′ ) = c 0 k 2 4 P(k -k ′ ), S(k) = |k|=|k ′ | s(k, k ′ )dk ′ (4.23)
where P(k) denotes P E 1 E 1 (k), i.e. the power spectral density for E 1 . As a matter of fact, in 1D media only two directions of propagation exist, thus k ′ = k or k ′ = -k. s(k, k ′ ) in eq. (4.23) reads:

s(k, k ′ ) =    c 0 k 2 4 P(0), k ′ = k c 0 k 2 4 P(2k), k ′ = -k , S(k) = c 0 k 2 4 P(0) + c 0 k 2 4 P(2k) (4.24)
Thus the scattering part in eq. (4.22) becomes:

|k|=|k ′ | s(k, k ′ )W ± (t, x; k ′ )dk ′ -S(k)W ± (t, x; k) = c 0 k 2 4 P(0)W ± (t, x; k) + c 0 k 2 4 P(2k)W ± (t, x; -k) - c 0 k 2 4 P(0) + c 0 k 2 4 P(2k) W ± (t, x; k) = c 0 k 2 4 P(2k)W ± (t, x; -k) - c 0 k 2 4 P(2k)W ± (t, x; k) (4.25)
where the two terms related to s(k, k) cancel each other out. Indeed, waves are scattered only in the direction opposite to their propagation direction.

Substituting eq. ( 4.25) into eq. (4.22) and using the convention defined by eq. (3.53), we get the 1D radiative transfer equations:

   ∂ t W + (t, x; -k) + c 0 ∂ x W + (t, x; -k) = c 0 k 2 4 P(2k) {W -(t, x; k) -W + (t, x; -k)} ∂ t W -(t, x; k) -c 0 ∂ x W -(t, x; k) = c 0 k 2 4 P(2k) {W + (t, x; -k) -W -(t, x; k)} (4.26)
Two new quantities emerge from adding and subtracting equations in eq. (4.26):

W s (t, x; k) = W + (t, x; -k) + W -(t, x; k), W d (t, x; k) = W + (t, x; -k) -W -(t, x; k) (4.27)
where W s and W d can be interpreted as the total energy density and the net forward-going density with respect to the propagation direction k |k| e x :

   ∂ t W s + c 0 ∂ x W d = 0 ∂ t W d + c 0 ∂ x W s = -c 0 k 2 2 P(2k)W d (4.28)
Equation (4.28) is equivalent to eq. (

Furthermore, by taking respectively the spatial and temporal derivatives of the second equation in eq. (4.28):

   ∂ 2 W d ∂ x∂t + c 0 ∂ 2 x W s = -c 0 k 2 2 P(2k)∂ x W d ∂ 2 t W d + c 0 ∂ 2 W s ∂t∂ x = -c 0 k 2 2 P(2k)∂ t W d (4.29)
and introducing the first equation of eq. (4.28) in eq. (4.29), we get:

   ∂ 2 t W s -c 2 0 ∂ 2 x W s = -c 0 k 2 2 P(2k)∂ t W s ∂ 2 t W d -c 2 0 ∂ 2 x W d = -c 0 k 2 2 P(2k)∂ t W d (4.30)
Two separate partial differential equations for W s and W d are obtained. If we consider a medium with no scattering, a same form as the traditional 1D wave equation is obtained in terms of W s or W d . For a given frequency k, eq. (4.30) can be seen as "damped" wave equations with a friction coefficient

c 2 0 k 2 2 P(2k) > 0.
It is known that the general solutions of the damped wave equation have an exponential decay over time with the coefficient e -Ct (C > 0 constant) [START_REF] Burq | Exponential decay for the damped wave equation in unbounded domains[END_REF][START_REF] Hancock | Solution to problems for the 1-d wave equation[END_REF][START_REF] Lin | The systems of second order partial differential equations with constant coefficients[END_REF]. The difference between wave energies in two opposite directions decreases quickly with time due to the conversion of energy in these two directions and the sum of energies in two directions decreases at the same rate as they spread out in space with the increase of time.

Scattering mean free path

In this section, an important characteristic length called the scattering mean free path l s for the radiative transfer regime is presented. The analytical solutions of the radiative transfer equation in 1D media are given firstly. Then we give the definition of l s and discuss its influence on this regime.

Let us assume two initial conditions for eq. (4.26) defined by

W + (t = 0, x; -k) = W + 0 (x; -k), W -(t = 0, x; k) = W - 0 (x; k) (4.31)
Equation (4.26) can be seen as the telegraph equation for a given wave number k. Thus its analytical solutions for a given k can be found in the following [START_REF] Hadeler | Reaction telegraph equations and random walk systems. Stochastic and spatial structures of dynamical systems[END_REF][START_REF] Howe | A kinetic equation for wave propagation in random media[END_REF][START_REF] Van Wijk | 1d energy transport in a strongly scattering laboratory model[END_REF]:

     W + (t, x; -k) = W + 0 (x -c 0 t)e -c 0 ls t + x+c 0 t x-c 0 t κ(t, x, y)W - 0 (y)dy + x+c 0 t x-c 0 t κ + (t, x, y)W + 0 (y)dy W -(t, x; k) = W - 0 (x -c 0 t)e -c 0 ls t + x+c 0 t x-c 0 t κ(t, x, y)W + 0 (y)dy + x+c 0 t x-c 0 t κ -(t, x, y)W - 0 (y)dy (4.32) with: l s (k) = 4 k 2 P(2k) , κ(t, x, y) = e -c 0 ls t 2l s I 0   c 2 0 t 2 -(y -x) 2 ) l s   , κ ± (t, x, y) = e -c 0 ls t 2l s I 1   c 2 0 t 2 -(y -x) 2 ) l s   c 0 t ∓ (y -x) c 2 0 t 2 -(y -x) 2 , (4.33)
where I 0 and I 1 are respectively the zero-order and the first-order modified Bessel functions of the first kind [START_REF] Furnas | Evaluation of the modified bessel function of the first kind and zeroth order[END_REF].

In fact, the first term of the two solutions in (4.32) represents the original coherent wave supplied by initial conditions, which decays exponentially with c 0 t over a characteristic length scale l s . Thus l s is the so called mean free path for scattering. It measures the distance of the exponential decrease of the energy density of coherent waves due to scatterings or successive interactions with the underlying medium. Roughly speaking, the second and the third terms in (4.32) with the modified Bessel functions describe the energy of incoherent waves which does not propagate ballistically. The second term results from the scattering of the original coherent waves in the opposite direction. The last term results form the conversion of coherent waves in the same direction when the coherent waves in the first term decrease.

When cumulative scatterings become significant during the propagation, the coherent energy becomes negligibly small. The incoherent energy dominates the field. W + and W -tend to a local equipartition between each other [START_REF] Howe | A kinetic equation for wave propagation in random media[END_REF]:

W + ≈ W -≈ e -x 2 2lsc 0 t 2 √ 2πl s c 0 t (4.34)
An example of analytical solutions for a single frequency wave is given here so as to illustrate the analysis above, with initial conditions defined by:

W + 0 (x; k) = 200δ (x)δ (k + k 0 ) m 3 s, W - 0 (x; k) = 100δ (x)δ (k -k 0 ) m 3 s (4.35)
where parameters are given by: k 0 = 157 m -1 , c 0 = 4 km/s. Thus l s = 3.46 m and the sinc 2 function is chosen as the power spectral density. 4.3d give their natural logarithmic values (denoted by ln) along x = ±c 0 t + ∆x in terms of c 0 t/l s (0 < ∆x ≪ 1). In fact, under the chosen initial conditions, it is known that W ± is a dirac function δ (0) on the characteristics x = ±c 0 t (eq. (4.32)), so we compute ln(W ± ) on a adjacent parallel line of the characteristics. The exponential decay of initial energies along the characteristic line (slope = ±c 0 ) is verified more rigorously. After the propagation over several l s , the coherent energy in x = ±c 0 t decreases quickly and it almost totally transforms into the incoherent energy centred around x = 0. After a long propagation time, W + and W -become more and more spread out in space. The distribution of total energy satisfies a local equipartition between all possible modes: in fig. 4.3e and fig. 4.3f, W + and W -become equal and symmetric about x = 0. Furthermore, note that l s depends on wave number k and the formula of power spectral density according to its definition (4.33). The example above is discussed for a given l s and a given frequency. Now we discuss its dependence on frequency. Choosing the three covariance functions: sinc 2 , exponential and gaussian presented in the table 4.1 as examples, their power spectral densities and scattering mean free paths are respectively computed and illustrated in fig. 4.4a and fig. 4.4b. Although their P(2k) and l s (k) have different supports in frequency, some common properties are found. For the high frequency limit, l s is approximatively inversely proportional to k 2 . For the low frequency limit, l s is approximatively inversely proportional to k 4 (for all P(k) sufficiently regular at k = 0 here). This phenomenon can be noticed in fig. 4.4b for the given parameters ε = 0.04, l c = 6 mm: l s tends to infinite for these two limits. Hence, no scattering of energy exists in either limit. Physically, this is due to the fact that in the low frequency domain (or the long wavelength limit) the medium becomes effectively homogeneous. In the high frequency limit, geometric optics becomes a good description for the classical wave propagation where the effect of scattering is expected to saturate. More studies of the influence of different correlation models on scattering mean free path can be found in [START_REF] Baydoun | Scattering mean free path in continuous complex media: Beyond the helmholtz equation[END_REF][START_REF] Khazaie | Influence of the statistical parameters of a random heterogeneous medium on elastic wave scattering: theoretical and numerical approaches[END_REF][START_REF] Van Der Baan | Acoustic wave propagation in one dimensional random media: the wave localization approach[END_REF].

Energy-based error estimation for HF wave in heterogeneous media Finally, note that what we discussed here remains valid when the waves are scattered by discrete scatterers that are randomly distributed in the medium [START_REF] Ryzhik | Transport equations for elastic and other waves in random media[END_REF][START_REF] Van Wijk | 1d energy transport in a strongly scattering laboratory model[END_REF].

Localization phenomenon in 1D heterogeneous media

Localization is a characteristic phenomenon when waves propagate in random media. It means that wave energy remains in a fixed bounded region of space close to the source at all times. It is well known that for (infinite) systems described by one spatial dimension, both quantum mechanical and classical waves are localized even if an infinitesimal amount of randomness is present for sufficiently large propagation distances [START_REF] Van Der Baan | Acoustic wave propagation in one dimensional random media: the wave localization approach[END_REF]. When waves are localized and the energy does not propagate, the radiative transfer regime is no more valid to describe the evolution of wave energy. Sheng [START_REF] Sheng | Introduction to wave scattering, localization and mesoscopic phenomena[END_REF] estimated that in 1D the localization length is approximately four times the mean free path (the factor four is not very accurate and valid especially for low frequency, but they nevertheless tell us that in 1D, the localization length is directly proportional to the mean free path). This offers the possibility of an intermediate range when the radiative transfer might hold.

In this section, we discuss a wave propagation problem in two heterogeneous media with different magnitudes of fluctuations. The localization phenomenon is then observed and analyzed.

Consider here the wave propagation in a heterogeneous bar with a constant density and a weakly fluctuating Young's modulus defined in eq. (4.21). A ricker force F(t) with a magnitude of 1 kN is applied in the centre of bar x = 0 so the wave front propagates in two opposite directions (fig. 4.5). Parameters are given in the table 4.2. The covariance function is chosen as sinc 2 with the parameters defined in the table 4.2. Due to the boundedness of its spectrum (fig. 4.4a), i.e. 2k/( 2π l c ) ≤ 1 for nonzero value, we take 2 × 3k max = 2π/l c (related to the cutoff frequency f c = 3 f max of ricker). Thus the power spectral density of the main frequency range of ricker does not equal zero, i.e. they have a relatively small and finite scattering mean free path (fig. 4.6). Consequently, l c = π 3k max = λ max 6 (it satisfies also the condition that λ is of the same order of magnitude as l c ). We choose two values of the variance σ 2 = ε for random media to observe the different phenomena. The mean free paths l s in these two cases are given in fig. 4.6. is considered here, as Bal and Pinaud [START_REF] Bal | Accuracy of transport models for waves in random media[END_REF] did in 2D random media. On the one hand, the existence of localization phenomena can be verified by observing the evolution of energy. On the other hand, it allows verifying the computation of Wigner transform as we introduced in the section 1.2.4.

ρ (kg/m 3 ) E 0 (GPa) σ 2 (E) = ε l c (mm) T R (ms) k max (m -1 ) λ max (m) h (
Firstly, W a is obtained analytically with eq. (4.31) subjected to some initial conditions. In this case, it should be noticed that we cannot impose analytically initial conditions for W a that is equivalent to the energy density of wave fields. However, we can choose an instant t 0 (after the source has vanished) and approximate the initial conditions W a (t = t 0 , x; k) by Wigner transform of numerical waves:

W a (t 0 , x; k) = R W ε [u h∆t ](t 0 , x; ω, k)dω (4.36)
Then the total energy densities (sum of kinetic and potential energies) in the time-space domain computed by analytical Wigner measure and numerical wave solutions in 1D are respectively defined as:

E[W a ](t, x) := 1 2 R ρω 2 W + a (t, x; k) + W - a (t, x; k)) + Ek 2 (W + a (t, x; k) + W - a (t, x; k) dk (4.37a) E[u h∆t ](t, x) := 1 2 R 2 {ρW ε [∂ t u h∆t ](t, x; ω, k) + EW ε [∂ x u h∆t ](t, x; ω, k)} dωdk (4.37b) = 1 2 ρ(∂ t u h∆t (t, x)) 2 + E(∂ x u h∆t (t, x)) 2
Equation (4.37a) is already introduced in the lemma 3. We can compute the evolution of the total energy with time, defined by:

Ẽ[W a ](t) = [-2,2] E[W a ](t, x)dx, E[ Ẽ[u h∆t ]](t) = E [-2,2] E[u h∆t ](t, x)dx (4.38)
so as to compare them and also to confirm the existence of localization. In fig. 4.10, it can be observed that the total energy Ẽ[W a ](t) is underestimated (red line) as we just explained. Thus we propose a "corrected" α Ẽ[W a ](t) with a correction factor:

α = E[ Ẽ[u h∆t ]](0) Ẽ[W a ](0) ≈ 1.2 (4.39)
that allows roughly eliminating the errors due to the estimated initial conditions. The comparison between the energy of averaged numerical waves (blue line) and the corrected energy of analytical Wigner measure (red dashed line) shows that wave energies are partly localized. T (s) According to the discussion above, the medium with the chosen amplitude of fluctuations ε = 0.04 leads to an obvious localization phenomenon during the whole observation period. Now we consider a smaller value ε = 0.008 and other parameters in the table 4.2 remain unchanged, then the same analysis is performed to verify the existence of localization and the validity of radiative transfer regime.

×10 -3 0 0.5 1 1.5 2 ×10 -7 E[ Ẽ[u]](t) Ẽ[W a ](t) α Ẽ[W a ](t)
Under this condition, the scattering mean free path is obtained in terms of frequency (fig. 4.6, l s (k max ) = 4.7 m). More precisely, l s is five times larger than that in the case with ε = 0.04. We still choose the same mesh size and time step as the example above for finite element solutions of wave equation. u h∆t for one realization is illustrated in fig. 4.13. With these two examples it is found that the smaller amplitude ε of fluctuations should be chosen to reduce the localization of energy so that the radiative transfer regime remains valid in a larger region. In the next section, the example with ε = 0.008 will be used for the computation of Wigner transform and residual errors.

Energy-based residual errors of numerical solutions of wave equation

As that defined in eq. (3.49) of residual errors for wave fields based on transport equation, here we can also define it with radiative transfer equations (4.11) in terms of W ε [u u h∆t ]. In this section, we only introduce this definition for 1D media for simplicity.

Theoretical residual errors based on 1D radiative transfer equation

   R + [u h∆t ] := ∂ t W + ε [u h∆t ] + c 0 ∂ x W + ε [u h∆t ] -c 0 l s (W - ε [u h∆t ] -W + ε [u h∆t ]) R -[u h∆t ] := ∂ t W - ε [u h∆t ] -c 0 ∂ x W - ε [u h∆t ] -c 0 l s (W + ε [u h∆t ] -W - ε [u h∆t ]) (4
     R sd T [u h∆t ] := ∂ t W s [u h∆t ] + c 0 ∂ x W d [u h∆t ] R sd T S [u h∆t ] := ∂ t W d [u h∆t ] + 2c 0 l s W d [u h∆t ] + c 0 ∂ x W s [u h∆t ] (4.41)
It can be noticed that R sd T contains only transport part of the radiative transfer equation and R sd T S contains both transport and scattering parts.

Either of these definitions as a representation of errors can be used since in fact we have R sd T = R + + R -and R sd T S = R + -R -. In latter numerical applications, R sd T and R sd T S are chosen to illustrate the results.

Numerical computation of Wigner transform of wave fields in random media

Before calculating numerically the defined errors in last section, the numerical computation of Wigner transform in random media is firstly specifically presented here.

Generally in random media, we evaluate physical quantities by taking the average over the ensemble of realizations of random media (as we did for the comparison of energy in the section 4.1.4). In addition, in the high frequency limit ε → 0, formal arguments show that the ensemble average of Wigner transform with respect to the randomness converge weakly to the Wigner measure, i.e. the solution of the radiative transfer equation in random media [START_REF] Bal | Self-averaging of wigner transforms in random media[END_REF][START_REF] Erdős | Linear boltzmann equation as the weak coupling limit of the random schrödinger equation[END_REF][START_REF] Spohn | Derivation of the transport equation for electrons moving through random impurities[END_REF]. Theoretically, W ε should be averaged over all possible realizations of the heterogeneous medium with given statistics, which causes high computational costs in numerical applications. In practice, we calculate, for example as in the last section, wave solutions in ten realizations of random media. A local average in the time-space domain is added within each realization so as to obtain an equivalent Energy-based error estimation for HF wave in heterogeneous media ensemble average, denoted always by E[W ε ], i.e.

E[W

ε ] = 1 j x N tx ∑ j ∑ (t k ,x k )∈I 1 ×Ω 1 ,k=1,..,N tx W ( j) ε (t k , x k ) (4.42)
Here I 1 and Ω 1 are respectively the time subinterval and subspace. j x is the total number of realizations and N tx is the total number of time-space points chosen in

I 1 × Ω 1 . W ( j)
ε is the Wigner transform of j-th realization of the random displacement field. j = 1, ..., 10 ( j x = 10) in our work. The chosen time-space points should be uncorrelated. In fig. 4.16, a simple illustration of a group of random points with a mean distance d m = 6l c between two adjacent points in t-x space. Then we can compute eq. ( 4.42) over these points and over ten realizations. 4.18a that energies are not only distributed in ω = ±c 0 k but also between these two characteristic lines. As we explained above, the latter are cross terms due to the interference of waves in two opposite directions. This phenomenon is more obvious than on fig. 4.17a since wave energies in two directions are comparable. Figure 4.18b gives E[W ε [u h∆t ]], an average in time-space domain and over 10 realizations. We find that energies are more concentrated along ω = ±c 0 k. Besides, it takes almost positive values over the entire phase space domain. Indeed, it can be understood by the energy properties of Wigner transform:

R×R W ε [u](t, x; ω, k)dtdx = | u(ω, k)| 2 > 0 (4.43) It is also observed that E[W + ε [u h∆t ]] and E[W - ε [u h∆t ]] are almost equal. In fact, E[W + ε [u h∆t ]], E[W - ε [u h∆t ]
] will always reach an equilibrium after a long time propagation. This result is explained in detail with the analytical solution of radiative transfer equation in the section 4.1.3.

In conclusion, both the ensemble average and the average over a set of uncorrelated space-time points are useful to reduce numerical fluctuations and cross terms in the computation of Wigner transform in random media. From now on, we still denote the average E[W ε [u h∆t ]] by W ε [u h∆t ] for simplicity.

Numerical application of residual errors based on radiative transfer equation

In this section, u h∆t with a relatively fine mesh size obtained in the numerical example (ε = 0.008) in the section 4.1.4 is used for calculating errors. Based on the calculations of Wigner transform for coherent waves and incoherent waves in the section 4.2.2, residual errors of u h∆t in eq. (4.41) are respectively computed and analyzed. In each case, we present firstly the errors R sd T and R sd T S in function of number of realizations N R . It is expected that errors are close to zero since a relatively fine mesh is used. Above all, in all figures that represent errors in this chapter, we do not consider the low frequencies since Wigner transform is computed in a bounded region and waves with a large wavelength cannot be totally and exactly captured. More precisely, we choose to observe only waves with:

k ≥ 2 × 2π ζ x = 40 m -1 (4.44)
Recall that ζ x is the computation region in space defined in the section 3.1.4 (fig. 3.9a) and we chose

ζ x = 8λ max in our work.
Firstly, we compute the residual errors for coherent forward waves in fig. 4.19. As we discussed in the section 4.2.2, a group of random points (N tx ≈ 300) in the time-space domain is chosen in the region around x = 0.5 m. R sd T and R sd T S are computed and averaged over these points and over ten random realizations. h∆t . And we find that it is relatively small compared to the total energy. It is observed high varying value of errors for lower frequency with different realizations (blue lines) and a local maximum of errors located at the maximum frequency k max = 157 m -1 for the averaged errors (red line).

Besides, we give here also |R + | and |R -| so as to analyze errrors for waves in two directions. Figure 4.21 gives absolute values of the residual errors in terms of |k|, normalized by the total energy. It is found that much smaller errors are obtained for |R + |. As we explained above, waves are mainly transmitted in the forward direction and only a small amount of waves are scattered to the backward direction, so errors for W - ε [u h∆t ] are polluted by numerical fluctuations of W + ε [u h∆t ] due to the frequency leakage and cross terms as we explained in fig. 4.17. 4.23, the related errors still have a main peak for low frequency. On the one hand, the cross terms cannot totally be suppressed by this average method. On the other hand, the numerical differentiation in eq. (4.41) for small and fluctuating energy density of low frequency will bring some additional errors (average may also reduce a part of such errors but very limited). These additional errors exist also for high frequencies but they are much smaller since we have much less cross terms. Furthermore, Wigner transform and the related residual errors can also be computed along the direction of propagation x = ±c 0 t. The advantage of this transformation of coordinates is that the edge effects or the discontinuities of boundary for numerical computation of Wigner transform may be smaller in this case. The numerical results can be found in Appendix F. It is found that a similar result is obtained and large errors are still found for low frequencies.

Filtering of the energy-based residual errors

Considering that the errors found in the last section are highly fluctuating and located largely at lower frequencies due to the superposition of plenty of cross terms in Wigner transform of a multi-frequency signal in random media, a filtering of residual errors is derived analytically in this section. Then we compute the new residual errors of the FE numerical solutions of wave equation in the last section 4.2.3.

We derive firstly the formula of filtered residual errors. In brief, it is based on the convolution of Wigner transform of wave fields with Wigner transform of a window function h(t, x).

where P (2) (r) is the second-order derivative of the covariance function P(r).

Proof. We can derive the 1D filtered radiative transfer equation by convolutions with eq. (4.26) for each term of the original radiative transfer equation.

For the two transport terms, we can prove the following result by lemma 3.4:

∂ t W ε [u ε ] * * W ε [h t h x ] = 2Re {W ε [∂ t u ε , u ε ]} * * W ε [h t h x ] (4.49)
As the auto-Wigner transform is always real and by the property 8, we have

∂ t W ε [u ε ] * * W ε [h t h x ] = 2Re {W ε [∂ t u ε , u ε ] * * W ε [h t h x ]} = 2Re {STFT[v ε , h t h x ]STFT * [u ε , h t h x ]} (4.50)
In the same way, it can be derived that:

∂ x W ε [u ε ] * * W ε [h t h x ] = 2Re {STFT[∂ x u ε , h t h x ]STFT * [u ε , h t h x ]} (4.51)
For the two scattering terms, we have:

k 2 4 P(2k)W ε [u ε ] * * W ε [h t h x ] = 1 4 1 (2π) 4 R 6 u ε t ′ - ετ 2 , x ′ - εy 2 u * ε t ′ + ετ 2 , x ′ + εy 2 h t t -t ′ - εs 2 h x x -x ′ - εz 2 h * t t -t ′ + εs 2 h * x x -x ′ + εz 2 R e -i(s-τ)ω ′ dω ′ :=(a) R k ′2 P(2k ′ )e -i(z-y)k ′ dk ′ :=(b) e i(sω+yk) dt ′ dx ′ dτdydsdz (4.52)
Firstly, the integration (a) can be seen as the Fourier transform of a function in the form f (x) = e iax (a constant). It is known that:

(a) = R e iτω ′ e -isω ′ dω ′ = 2πδ (s -τ) (4.53)
Secondly, for the integration (b), by using respectively the definition of the inverse Fourier transform and the Fourier transform of function f

(x) = x n [35]: R k n e -jxk dk = 2πi n δ (n) (x) (4.54)
where n is a natural number and δ (n) (x) is the n-th distribution derivative of the dirac delta function, we can derive:

Energy-based error estimation for HF wave in heterogeneous media

(b) = R k ′ 2 R P(x)e -i2xk ′ dx e -i(z-y)k ′ dk ′ = R P(x) R k ′ 2 e -i(2x+z-y)k ′ dk ′ dx = R P(x) 2πi 2 δ (2) (2x + z -y) dx (4.55)
Then using the derivative property of the dirac delta function [START_REF] Ghatak | The dirac delta function[END_REF]:

f (x)δ (n) (x)dx = - ∂ f ∂ x δ (n-1) (x)dx (4.56)
and we get:

(b) = -2π R P(x)δ (2) (2x + z -y) dx = - π 4 P (2) y -z 2 (4.57)
Introducing eq. (4.53) and eq. (4.57) to eq. ( 4.52), we get:

k 2 4 P(2k)W ε [u ε ] * * W ε [h t h x ] = 1 4 1 (2π) 4 R 6 u ε t ′ - ετ 2 , x ′ - εy 2 u * ε t ′ + ετ 2 , x ′ + εy 2 h t t -t ′ - εs 2 h x x -x ′ - εz 2 h * t t -t ′ + εs 2 h * x x -x ′ + εz 2 2πδ (s -τ) (- π 4 )P (2) y -z 2 e i(sω+zk) dt ′ dx ′ dτdydsdz = - π 16
1 (2π) 3 R 5 u ε t ′ - ετ 2 , x ′ - εy 2 u * ε t ′ + ετ 2 , x ′ + εy 2 h t t -t ′ - ετ 2 h x x -x ′ - εz 2 h * t t -t ′ + ετ 2 h * x x -x ′ + εz 2 P (2) y -z 2
e i(τω+zk) dt ′ dx ′ dτdydz (4.58) Then changing the variables in the equation above (t ′ , τ) by (α, β ) and (x ′ , y, z) by (m, n, l) with:

t -t ′ + ετ 2 = εα t -t ′ -ετ 2 = εβ ,        x -x ′ + εz 2 = εm x -x ′ -εz 2 = εn y-z 2 = l (4.59)
thus R 2 dt ′ dτ = εdαdβ and R 2 dx ′ dydz = 2εdmdndl. We have: Same derivation can be applied for the scattering term in eq. (4.26) with "-k" (energy in the opposite direction):

k 2 4 P(2k)W ε [u ε ] * * W ε [h t h x ] = - ε 2 16 1 (2π) 2 R R 2 u ε (t -εα, x -εl -εm) h * t (εα)h * x (εm) e i(αω+mk) dαdm R 2 u * ε (t -εβ , x -εl + εn) h t (εβ )h x (εn) e -i(β ω+nk) dβ dn P (2) (l)dl = - 1 16 R STFT[u ε , h t h x ](t, x -εl; ω, k)STFT * [u ε , h t h x ](t, x + εl; ω, k) P (2) (l)dl
k 2 4 P(2k)W ε [u ε ](-k) * * W ε [h t h x ] = - ε 2 16 1 (2π) 2 R R 2 u ε (t -εα, x -εl -εm) h * t (εα) h x (εm) e i(αω-mk) dαdm R 2 u * ε (t -εβ , x -εl + εn) h t (εβ ) h * x (εn) e -i(β ω-nk) dβ dn P (2) (l)dl = - 1 16 R STFT[u ε , h t h * x ](t, x -εl; ω, -k)STFT * [u ε , h t h * x ](t, x + εl; ω, -k) P (2) (l)dl (4.61)
Combining the two transport terms in equations 4.50, (4.51) and the two scattering terms in eqs. (4.60), (4.61) and projecting in "±", the filtered radiative transfer equation for forward and backward waves is obtained:

2Re {STFT[v h∆t , h t h x ](t, x; c 0 k, ∓k)STFT * [u h∆t , h t h x ](t, x; c 0 k, ∓k)} ± 2c 0 Re {STFT[∂ x u h∆t , h t h x ](t, x; c 0 k, ∓k)STFT * [u h∆t , h t h x ](t, x; c 0 k, ∓k)} = - 1 16 R STFT[u h∆t , h t h * x ](t, x -εl; c 0 k, ±k)STFT * [u h∆t , h t h * x ](t, x + εl; c 0 k, ±k) P (2) (l)dl + 1 16 R STFT[u h∆t , h t h x ](t, x -εl; c 0 k, ∓k)STFT * [u h∆t , h t h x ](t, x + εl; c 0 k, ∓k) P (2) (l)dl (4.62)
Then the filtered residual errors are obtained in eq. (4.48) as we did in eq. (4.41).

* * W[h t h x ] , (b) R sd T S * * W[h t h x ]
We consider the same wave propagation problem in ten realizations of random media as in section 4.2.3. Then we compute the absolute value of the filtered residual errors R sd

T [u h∆t ] * * W ε [h t h x ] and R sd T S [u h∆t ] * * W ε [h t h x ]

Comparison of residual errors for two refinements

In this section, in order to study the behaviour of the proposed residual errors with respect to the discretization size in space and in time, numerical wave solutions obtained by using a space mesh twice coarser and a time step twice larger than those used in the section 4.2.3 are considered, so we have λ (3k max )/h = 10, ∆t = h/c 0 . Otherwise, all the numerical parameters for calculating the energy-based residual errors, such as the computation region of Wigner transform, the random points chosen for average, remain unchanged. The residual errors calculated for both discretization sizes are compared and it is hoped that they would decrease when the discretization in space and in time is refined.

The residual errors defined in eq. (4.41) and the filtered residual errors defined in eq. (4.48) are computed respectively for coherent waves and incoherent waves as we did in the section 4. ), as expected both residual errors and filtered residual errors are smaller, especially over the main frequency range defined by the considered ricker signal, for the wave solutions calculated with the finer discretization. Nevertheless, the decrease of residual errors between the two discretizations is largely reduced after the filtering, which is not yet well understood for the moment. As far as the incoherent waves are concerned, both residual errors and filtered residual errors give unsatisfying results, in fact no obvious decrease of errors is observed with the finer discretization (figs. 4.28 and 4.29). There are several reasons that may account for this phenomena. In section 3.3.1, it was observed that numerical dispersion and dissipation may influence the quality of estimated errors. Here after a long propagation time, the dispersion or dissipation of incoherent waves become significant. Thus the decrease of errors are not found as we observed in fig. 3.13b. Besides, when we refine the mesh size and the time step, the elastic parameter E is also refined. It means that the two solutions here are in fact calculated with two numerical models that have slightly different physical properties. That may also disturb the results of errors. There exist also other possible reasons, such as numerical fluctuations resulting from frequency leakage, numerical differentiation and the distribution of cross terms, existence of localization, etc. Further studies are necessary to understand these results and to resolve this problem.

Conclusion

This chapter dealt with energy-based residual error for high frequency wave fields in heterogeneous media. We introduced the heterogeneous media modelled with realizations of random fields of mechanical properties. Radiative transfer equation that describes the transport and the scattering of wave energy is obtained in these media in the weak coupling regime. 1D radiative transfer equations were studied and solved analytically. The characteristic length for wave propagation in random media, i.e. the scattering mean free path was introduced. It was observed that in low frequency and high frequency limits, the scattering mean free path becomes unbounded. The localization phenomenon that always exists in 1D was analyzed by considering the wave propagation problems in two random media with different amplitudes of fluctuations. This analysis allowed finding a region in a random medium where the radiative transfer regime still holds. Then based on the radiative transfer equation, the residual errors of numerical solutions of wave equation were defined and computed numerically. It was found that the obtained errors are very small compared to the total energy. Also, large part of errors are concentrated around low frequencies due to cross terms of Wigner transform of random fields.

Conclusion 105

In order to eliminate this effect, a filtering of residual errors was derived based on the convolution properties of Wigner transform. It also favours the smoothing of the results. Numerical results showed that this filtering with window function can reduce partly the pollution on errors caused by cross terms, especially for low frequency. Finally, the decrease of errors with refinement was found with coherent waves as we expected. However, the influence of filtering and the unexpected results for incoherent waves remain to be improved.

Chapter 5

Conclusions and perspectives

Wave propagation in heterogeneous media, especially in the high frequency regime is always a challenging and complex problem. Its numerical modelling and simulation require very accurate numerical methods with high computational and memory costs. This work is dedicated to the development of tools of error estimation for finite element solutions of wave equation in heterogeneous media, which will be used to find solutions with a desired precision at a moderate cost.

Firstly, an explicit a posteriori upper bound for the errors between the exact solution and the reconstructed approximate solution in an L ∞ -norm of energy in time was derived analytically for the elastodynamic equation.

The basic idea of this derivation is based on the residual error method and the use of reconstructions. In consideration of the requirement of regularities for the error bound, we have proposed two reconstructions in time with different quantities or orders for displacement and velocity fields, and two reconstructions in space for stress field: one is a basic reconstruction based on an element-wise weak equilibrium relation, another is the statically admissible reconstruction, inspired by Ladevèze [START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF]. We studied firstly the behaviour of the exact numerical errors in homogeneous and heterogeneous cases with two different degrees of heterogeneity. The finite element solutions are solved with different mesh sizes and time steps. We found that the errors are decreasing with the refinement in all cases. An exchange between errors in kinetic and potential energy norms has been observed when waves interact with boundaries and interfaces. Finally, by choosing different combinations of reconstructions, the estimated errors were compared with the exact errors in the homogeneous case. A fully computable error upper bound has been found with our error estimator but the ratio of the estimated and exact errors is increasing as the computational effort grows. It means that the asymptotic accuracy of the obtained estimator remains to be improved. Some discussions and analysis on the results showed that the residual errors related to the time derivative of the equilibrium equation are dominant in the contribution of total estimated errors.

Secondly, considering the difficulties encountered in numerical applications of the error bound for the elastodynamic equation, a specific case, the high frequency wave propagation in the weak coupling Conclusions and perspectives regime was studied. As the traditional error estimators on the wave equation may not be efficient in this regime, we proposed a quantification of errors with a kinetic model in terms of energy densities, i.e. the transport equation in homogeneous media and the radiative transfer equation in heterogeneous media. This kinetic model allows smoothing out the highly oscillatory features of waves. Based on this new equation, we defined analytically a new residual-type error estimation method. The errors are quantified in terms of energy quantities of numerical waves, which are calculated by the Wigner transform of wave fields.

In 1D homogeneous media, the localization of errors around the characteristics in phase space was observed. Considering the existence of extra fluctuations by frequency leakage, a filtering of errors with window functions in phase space was proposed and a global residual error in time-space domain was defined. The decrease of local and global errors when the discretization size decreases was observed, which allows validating our definition of errors and methods of calculation.

In 1D heterogeneous media, the residual errors in terms of the total energy and the net forwardgoing energy were proposed analytically. In numerical applications, a 1D medium with a small magnitude of fluctuations was firstly found so that less localization exists and the radiative transfer regime holds. Also an agreement between the radiative transfer energy and the numerical wave energy allowed validating our calculations. Then the defined errors in the chosen medium with a fine mesh were computed with the energy quantities of numerical wave solutions. The results revealed that the defined errors are relatively small compared to the total energy, that allowed validating our definition of errors. In order to reduce the influence of cross terms and other fluctuations on the distribution of errors in terms of frequency and to obtain a smoother result of errors, a filtering of residual errors in terms of the STFT of wave fields was derived. It is based on the convolution of radiative transfer equation with Wigner transform of a window function. This filtering allowed reducing especially the cross terms accumulated in low frequencies. Finally, the comparison of errors in phase space for numerical solutions with two different discretizations sizes was obtained. For coherent waves, residual errors are smaller for the finer discretization size as we expected. However, the decrease of residual errors between the two discretizations was reduced after the filtering. For incoherent waves, no obvious decrease of errors was found with the finer discretization. Some possible reasons were given and a detailed analysis remains to be made in future works. Some improvements, perspectives and difficulties are discussed here.

In a short term, efforts should be made firstly for resolving the existing problems. For the explicit error bound defined for elastodynamics that has a weak asymptotic accuracy, the reconstructions in time and in space should be reconsidered. They were defined with some desired regularities but they also bring in additional errors since they can be seen as a post-treatment of approximate solutions. Further efforts are needed to propose some new reconstructions so as to achieve a better equilibrium, especially for the time derivative of equilibrium on the entire time interval. For the energy-based residual errors in heterogeneous media, future works should be firstly focused on discussing and analyzing the phenomena in the comparison of errors with two discretization sizes as we discussed above. Besides, at present the identification of energy-based errors in terms of frequency seems difficult to achieve due to the superposition of cross terms between all frequency components. Since the filtering with the STFT allows already reducing largely cross terms in low frequencies, this idea may be applied in future works. It may allow quantifying the effect of numerical dispersion. In the whole work, we merely defined errors to quantify the difference between the exact solutions and the approximate solutions solved by time discontinuous Galerkin method. Other solvers or discretization methods, such as Newmark, should also be tested to observe their influence on the estimator. Note that the proposed reconstructions may be adapted to the choice of solvers. Furthermore, in the radiative transfer regime, only wave propagation in an infinite region is considered. We ignore the calculation of Wigner transform in the boundaries in the homogeneous case when waves are totally reflected by the boundary. Difficulties arise when the energy densities reach the boundaries of the domain or an interface. Some researchers have studied systematically the boundary or interface conditions for the high frequency waves propagation in media with some specific boundaries or interfaces [START_REF] Bal | Transport theory for acoustic waves with reflection and transmission at interfaces[END_REF][START_REF] Savin | Transient vibrational power flows in slender random structures: Theoretical modeling and numerical simulations[END_REF][START_REF] Staudacher | Conservative numerical schemes for high-frequency wave propagation in heterogeneous media[END_REF]. But generally it is difficult to construct the Wigner transform of wave fields with respect to boundary conditions in the case with arbitrary geometries of the boundary and interfaces. Future works may be focused on errors based on the initial boundary value problem for the radiative transport equation in some simple cases with specific boundaries and interfaces.

Finally, the energy-based residual errors were analytically defined in 3D, so future works may extend the numerical results of 1D media to a higher dimension. It should be noticed that localization exists always in 1D media so applications on higher dimension are necessary. Note that the computation and storage costs will be very expensive in higher dimensions in the sense that the Wigner transform depends on (t, x; ω, k), i.e. six independent variables in 2D, and eight independent variables in 3D. However, considering that the dispersion relation ω = c p,s |k| is always obtained and involved in the radiative transfer equation, we can study only the quantities on the hypersurface defined by this relation (2D media), that may allow reducing significantly the storage cost. However, for numerical solutions, we can never remove these fluctuations with either of these two methods because of the evolution in time of the numerical dispersions. Also, it cannot be applied in the heterogeneous case since the displacement is no more constant along the propagation direction.

Spectral leakage in FFT

W ε [u ε ](t, x; ω, k) = 1 (2π) 2 R 2 u ε 1 √ 2c 0 (χ ⊥ + χ ∥ ) - ε 2 1 √ 2c 0 (φ + ψ), 1 √ 2 (-χ ⊥ + χ ∥ ) - ε 2 1 √ 2 (-φ + ψ) u * ε 1 √ 2c 0 (χ ⊥ + χ ∥ ) + ε 2 1 √ 2c 0 (φ + ψ), 1 √ 2 (-χ ⊥ + χ ∥ ) + ε 2 1 √ 2 (-φ + ψ) e i 1 √ 2c 0 (φ +ψ)ω+ 1 √ 2 (-φ +ψ)k 1 c 0 dφ dψ = 1 (2π) 2 R 2 u ε 1 √ 2c 0 ((χ ⊥ - εφ 2 ) + (χ ∥ - εψ 2 )), 1 √ 2 (-(χ ⊥ - εφ 2 ) + (χ ∥ - εψ 2
))

u * ε 1 √ 2c 0 ((χ ⊥ + εφ 2 ) + (χ ∥ + εψ 2 )), 1 √ 2 (-(χ ⊥ + εφ 2 ) + (χ ∥ + εψ 2
))

ũ(χ ⊥ , χ ∥ ) := u 1 √ 2c 0 (χ ⊥ + χ ∥ ), 1 √ 2 (-χ ⊥ + χ ∥ ) ; Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) := 1 (2π) 2 R 2 ũε (χ ⊥ - εφ 2 , χ ∥ - εψ 2 ) ũ * ε (χ ⊥ + εφ 2 , χ ∥ + εψ 2
)e i(φ k χ ⊥ +ψk ∥ ) dφ dψ with the change of variables in the phase space

   k ⊥ = 1 √ 2 ( ω c 0 -k) k ∥ = 1 √ 2 ( ω c 0 + k)
. Then the projections of Wigner transform on the forward and backward directions can be expressed as:

W + ε = W ε [u ε ](t, x; ω = c 0 k, -k) = 1 c 0 Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ = √ 2k, k ∥ = 0) = 1 c 0 W+ ε W - ε = W ε [u ε ](t, x; ω = c 0 k, k) = 1 c 0 Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ = 0, k ∥ = √ 2k) = 1 c 0 W- ε (F.7)
The local wave-characteristic basis and the associated coordinate system are illustrated in fig.

∂ t W ε [u ε ](t, x; ω, k) = 1 c 0 ∂ χ ⊥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ )∂ t χ ⊥ + 1 c 0 ∂ χ ∥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ )∂ t χ ∥ = 1 c 0 ∂ χ ⊥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) c 0 √ 2 + 1 c 0 ∂ χ ∥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) c 0 √ 2 = 1 √ 2 ∂ χ ⊥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) + 1 √ 2 ∂ χ ∥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) (F.8)
and

∂ x W ε [u ε ](t, x; ω, k) = 1 c 0 ∂ χ ⊥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ )∂ x χ ⊥ + 1 c 0 ∂ χ ∥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ )∂ x χ ∥ = - 1 c 0 ∂ χ ⊥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) 1 √ 2 + 1 c 0 ∂ χ ∥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) 1 √ 2 = - 1 √ 2c 0 ∂ χ ⊥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) + 1 √ 2c 0 ∂ χ ∥ Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) (F.9)
Combing the above two equations and using eq. (F.6) and eq. (F.7), the two radiative transfer equations for W+ and Win eq. (4.26) become:

   √ 2∂ χ ∥ W+ [ ũε ] = 1 8 k 2 ⊥ P( √ 2k ⊥ ) W-[ ũε ] -W+ [ ũε ] √ 2∂ χ ⊥ W-[ ũε ] = 1 8 k 2 ∥ P( √ 2k ∥ ) W+ [ ũε ] -W-[ ũε ]
(F.10)

Theorem F.1. The residual errors in the local basis associated to the wave propagation characteristics are defined with: The filtering of residual errors can also be transformed under (χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) in the same way as we did in theorem 4.6. The detailed derivation is neglected and the results are given here. Then we can calculate eq. (F.12) and eq. (F.14) numerically for coherent waves and incoherent waves as we did in section 4.2.3.

   R+ [ ũε ] = √ 2∂ χ ∥ W+ [ ũε ] -1 8 k 2 ⊥ P( √ 2k ⊥ ) W-[ ũε ] -W+ [ ũε ] R-[ ũε ] = √ 2∂ χ ⊥ W-[ ũε ] -1 8 k 2 ∥ P( √ 2k

Fig. 1 .Fig. 1 .

 11 Fig. 1.1 (a) polycrystalline microstructure, (b) simulation of quasi-P wave propagation in polycristals (by MSSMat, CentraleSupélec)

Fig. 1 . 3

 13 Fig. 1.3 Body waves (left) and surface waves (right) (Encyclopaedia Britannica, Inc.)

Fig. 2 .

 2 Fig. 2.1 (a) finite element solutions (u n h∆t , v n h∆t ) in time, (b) exact and finite element stress in space

Fig. 2 . 3 Fig. 2 . 4

 2324 Fig. 2.3 Truncated ricker signal and its frequency spectrum (amplitude normalized to 1): (a) F(t), (b) |FFT (F)| (fast Fourier transform)

Fig. 2 . 6 1 2

 261 Fig. 2.6 Numerical results of velocity fields v(t, x) and local errors 1 2 ρe 2 v + E(∂ x e u ) 2 (t, x) (normalized by total energy L 0

  2.8b).

Fig. 2 . 7

 27 Fig. 2.7 Homogeneous media (λ /h = 32): (a) velocity fields v(t, x), (b) ∥e v ∥ M (t) and ∥e u ∥ K(Ω) (t)

Fig. 2 . 8

 28 Fig. 2.8 Heterogeneous media with one inclusion (λ 2 /h = 32): (a) velocity fields v(t, x), (b) ∥e v ∥ M (t) and ∥e u ∥ K(Ω) (t)

Fig. 2 . 9

 29 Fig. 2.9 Heterogeneous media with five inclusions (λ 2 /h = 32): (a) velocity fields v(t, x), (b) ∥e v ∥ M (t) and ∥e u ∥ K(Ω) (t)

Fig. 2 .

 2 Fig. 2.10 A log-log plot of exact errors ∥ev ∥ 2 M(Ω),∞(T f ) + ∥e u ∥ 2 K(Ω),∞(T f ) (black line is used as reference with slope -1)

Fig. 2 .

 2 Fig. 2.11 Error estimation in the homogeneous case: (a) exact errors and estimated errors (log-log graph), (b) effectivity index (semi-log graph)

Chapter 3 Energy

 3 -based a posteriori error estimation for high frequency wave equation in homogeneous media Contents 3.1 Wigner transform applied to the elastic wave equation . . . . . . . . . . . . . 43 3.1.1 Temporal Wigner transform in signal processing . . . . . . . . . . . . . 43 3.1.2 Discussion and analysis of the interference property of Wigner transform 45 3.1.3 Scaled spatio-temporal Wigner transform in the high frequency limit . . . 49 3.1.4 Discrete Wigner transform and analysis of parameters . . . . . . . . . . 51 3.2 Transport energy-based residual errors of numerical solutions of wave equation 58

  3) Energy-based error estimation for HF wave in homogeneous media Property 2 (Conjugation property). W[F, G] = (W[G, F]) * (3.4) Thus W[F] is real for all (t; ω). Note that in the case of Wigner transform of two vectorial fields F and G, W[F, G] = (W[G, F]) * T . Property 3 (Correlation of a sum property). The Wigner transform of a sum of two signals is the square of the correlation in time-frequency of their Wigner transforms.

Fig. 3 .

 3 Fig. 3.1 (a) real part of a complex-valued signal F(t) (eq. (3.11), eq. (3.15)) composed of two gaussian functions localized in both time and frequency domains, (b) Wigner transform in time-frequency domain W[F](t; ω): two auto-terms (positive, red round zone); one cross term (positive or negative, oscillating zone)

3 Fig. 3 . 2

 332 Fig. 3.2 Energy content of F(t) (composed of two disjoint gaussian functions): R W[F](t; ω)dω normalized to total energy

Fig. 3 . 3

 33 Fig. 3.3 Wigner transform of a signal G(t) (eq. (3.16), eq. (3.15)) composed of two real-valued gaussian functions in time-frequency domain W[G](t; ω): four auto-terms (positive, red round zone); six cross terms (positive or negative, oscillating zone)

  ζ t = N s ∆t is the sampling length. The summation variable τ in Wigner transform is discretized by τ/2 = M∆t, M ∈ Z. The discretization of frequency is defined as ω = 2π L N s ∆t , L ∈ Z. Then the discrete Wigner transform W[F][N; L] is defined as (we use the same notation W[F] for continuous and discrete Wigner transform):

Fig. 3 . 4

 34 Fig. 3.4 Ricker signal and its frequency spectrum: (a) F(t), (b) |FFT (F)|

Figure 3 .

 3 Figure3.5a -3.5d give the numerical results of W[F] in terms of time-frequency for four different values of ζ t . We find also that the maximum value of Wigner transform lies at zero frequency, as a result of superposition of cross terms between all ω 0 and -ω 0 (ω 0 ∈ [0, 2π f max ]) since ricker is a real-valued signal in time. The second maximum value of Wigner transform is located around the maximum frequency.In fig.3.5e, we compare the four plots fig.3.5a -3.5d for a fixed frequency f = f max . According to the convolution property of Wigner transform in frequency (eq. (3.9a)), the change of computation region, i.e. the size of rectangular window function has almost no influence on time resolution. The last three lines almost coincide, but ζ t = T max = T R /2 differs most from others because it contains only the half signal in time and has a low-resolution in time. It changes principally the frequency resolution. Figure3.5f shows their values for a fixed time instant t = 0: we observe that ζ t = 2T max already allows identifying a maximum frequency around f max = 500 Hz as expected. With increasing ζ t , a spectrum with higher resolution is obtained.

Fig. 3 . 5 Fig. 3 . 6 Fig. 3 . 7 Table 3 . 1

 35363731 Fig. 3.5 Wigner transform of ricker signal W[F](t; f ) for (a) ζ t = T max , (b) ζ t = 2T max , (c) ζ t = 4T max , (d) ζ t = 8T max ; comparison of (e) W[F](t; f = f max ) and (f) W[F](t = 0; f ) for the four considered values of ζ t

Fig. 3 .

 3 Fig. 3.8 W[u](t, x; ω, k) for (t 0 , x 0 ) at wave front with: (a) (ζ t , ζ x ) = (2T max , 2λ max ), (b) (ζ t , ζ x ) = (4T max , 4λ max ), (c) (8T max , 8λ max ); (d) their value at the line perpendicular to ω = -c 0 k for the maximum frequency k = -k max

Fig. 3 . 9 Table 3 . 2 Fig. 3 .Figure 3 .

 393233 Fig.3.9 Wave displacement: (a) u(t, x) in the whole computed time-space domain (with two computation regions for Wigner transform at two points A (t A , x A ) and B (t B , x B )), (b) u(t, x 0 ) for x 0 = 0.4 m, (c) u(t 0 , x) for t 0 = 0.1 ms

Figure 3 .

 3 Figure 3.11 give local residual errors distributed in phase space for the points A and B with different numerical solutions. It can be found the tendency of decreasing of errors with increasing refinement. Errors at point B are generally larger that at point A for same size of refinement due to the increase of errors with time.

Fig. 3 .Fig. 3 .Fig. 3 .

 333 Fig. 3.11 (a) R ± for point A in phase space for different refinements; (b) R ± for point B in phase space for different refinements

√ 2ζ x = 2 √

 2 2∆k. For backward waves, errors are distributed in the first quadrant, thus k ⊥ ∈ [0, 3 √ 2k max ].

Figure 3 .

 3 Figure 3.14 illustrates the effect of window function on residual errors for forward waves (for λ /h = 4 in fig. 3.11a).

Fig. 3 .Fig. 3 .

 33 Fig. 3.14 Filtering for residual errors in phase space

Fig. 3 .Fig. 3 .

 33 Fig. 3.16 Comparison of residual errors between numerical solutions and exact solutions: (a) R(t, x 0 ) for x 0 = 0.4 m, (b) R(t 0 , x) for t 0 = 0.1 ms

Contents 4 . 1

 41 Radiative transfer equation in heterogeneous media in the weak coupling limit 73 4.1.1 Statistical description of weak randomly-fluctuating media . . . . . . . . 73 4.1.2 Radiative transfer equation . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.1.3 Scattering mean free path . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.1.4 Localization phenomenon in 1D heterogeneous media . . . . . . . . . . 83 4.2 Energy-based residual errors of numerical solutions of wave equation . . . . 90 4.2.1 Theoretical residual errors based on 1D radiative transfer equation . . . . 90 4.2.2 Numerical computation of Wigner transform of wave fields in random media 91 4.2.3 Numerical application of residual errors based on radiative transfer equation 94 4.3 Filtering of the energy-based residual errors . . . . . . . . . . . . . . . . . . 97 4.4 Comparison of residual errors for two refinements . . . . . . . . . . . . . . 102 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

l 2 cFig. 4 .

 24 Fig. 4.1 (a) normalized sinc 2 , exponential and gaussian covariance functions, (b) sinc 2 function with different correlation lengths

Fig. 4 . 2

 42 Fig. 4.2 One realization of 2D random media generated by the covariance function sinc 2 with different correlation lengths: (a) l c = l 0 , (b) l c = 2l 0

Figure 4 .

 4 Figure 4.3a and fig. 4.3b give results of W ± in the time-space domain. It can be observed the exponential decay over propagation time-distance along x = ±c 0 t (red line). More precisely, fig. 4.3c and fig.4.3d give their natural logarithmic values (denoted by ln) along x = ±c 0 t + ∆x in terms of c 0 t/l s (0 < ∆x ≪ 1). In fact, under the chosen initial conditions, it is known that W ± is a dirac function

Fig. 4 . 3

 43 Fig. 4.3 Analytical solutions of W ± in the time-space domain: (a) W + (t, x) (red line: x = c 0 t), (b) W -(t, x) (red line: x = -c 0 t), (c) ln(W + ) on x = c 0 t + ∆x, (d) ln(W -) on x = -c 0 t + ∆x, (e) W + (•, x) for four time instants, (f) W -(•, x) for four time instants t = 0.4 ms, 2.6 ms, 4.4 ms and 7 ms, which correspond respectively to c 0 t/l s (k 0 ) = 0.3, 3, 5 and 8.

Fig. 4 . 4

 44 Fig. 4.4 Normalized power spectral density and mean free path for different correlation models: (a) P(2k), (b) l s (k) (parameters: ε = 0.04, l c = 6 mm)

Fig. 4 . 5

 45 Fig. 4.5 Model in 1D heterogeneous media

Fig. 4 . 6

 46 Fig. 4.6 Mean free path l s (k) (parameters: ε = 0.04, 0.008, l c = 6 mm)

Fig. 4 . 7

 47 Fig. 4.7 Wave displacement in time-(half) space for one realization (ε = 0.04)

6 .

 6 Fig. 4.8 Energy density in the time-space domain (t ∈ [0, 2 ms], x ∈ [-2 m, 2 m]) obtained by one realization of numerical wave fields E[u h∆t ]

Fig. 4 . 9

 49 Fig. 4.9 Energy density in the time-space domain (t ∈ [0, 2 ms], x ∈ [-2 m, 2 m]) obtained by (a) (b): analytical Wigner measure E[W a ], (c) (d): average of 10 numerical wave fields E[E[u h∆t ]]

Fig. 4 .Fig. 4 .

 44 Fig. 4.10 Comparison of the evolution of energy in time obtained by analytical/numerical solutions: E[ Ẽ[u h∆t ]] (blue solid line), Ẽ[W a ] (red solid line), α Ẽ[W a ] (red dashed line)

25 EFig. 4 .

 254 Fig. 4.12 Comparison of variation of energy in time obtained by analytical/numerical solutions for three instants: t 1 = 0.5 ms, t 2 = 0.75 ms, t 3 = 1.15 ms (final time = 2 ms)

Fig. 4 .

 4 Fig. 4.13 Wave displacement in time-(half) space for one realization (ε = 0.008)

  Fig. 4.14 Comparison of variation of energy in time obtained by analytical/numerical solutions for three instants: t 1 = 1 ms, t 2 = 3 ms, t 3 = 6 ms (final time = 8 ms)

Fig. 4 .

 4 Fig. 4.15 Comparison of total energy obtained by analytical/averaged numerical solutions in time for x = 0 (white dashed line: smoothed value of E[E(u h∆t )](t, 0)): (a) ε = 0.04, (b) ε = 0.008

Definition 4 . 2 .

 42 Based on 1D radiative transfer equation (4.26 and 4.28), residual errors in terms of Wigner transform of wave fields are defined in the following: Residual errors based on 1D radiative transfer equation in terms of W ± are defined as:

1 Fig. 4 .

 14 Fig. 4.16 Randomly chosen time-space points for averaging the Wigner transform of numerical waves

Fig. 4 .

 4 Fig. 4.17 Computation of Wigner transform for coherent forward waves in random media: (a) (b) W ε [u h∆t ] for one point (t 0 , x 0 ) with two different scales, (c) (d) E[W ε [u h∆t ]] averaged in time-space and in ten realizations with two different scales

Figure 4 .Fig. 4 .Figure 4 .

 444 Figure 4.17 gives numerical results of W ε [u h∆t ](t, x; ω, k) of coherent forward wave fields. Similar to the homogeneous case in fig. 3.10a, energies are mainly distributed along ω = -c 0 k since the coherent energy are dominant. However, some fluctuations are also observed here (fig. 4.17b) compared to fig. 3.10a. On the one hand, more discontinuities in the boundary of computation region

  zoom

Fig. 4 .

 4 Fig. 4.19 Computation area presented with u(t, x)

Figure 4 .

 4 Figure 4.20 gives absolute values of the residual errors R sd T and R sd T S in terms of |k| for three random realizations and the average of ten realizations, normalized by the total energy E u 2h∆t . And we find that it is relatively small compared to the total energy. It is observed high varying value of errors for lower frequency with different realizations (blue lines) and a local maximum of errors located at the maximum frequency k max = 157 m -1 for the averaged errors (red line).

Fig. 4 .Fig. 4 .Fig. 4 . 4 . 3 97 Figure 4 .

 44443974 Fig. 4.20 Residual errors along ω = ±c 0 k for coherent wave: (a) R sd T , (b) R sd T S

Fig. 4 .Figure 4 .

 44 Fig. 4.25 Residual errors along ω = ±c 0 k for incoherent wave: (a) |R + |, (b) |R -|

(4. 60 ) 4 . 3

 6043 Filtering of the energy-based residual errors 101

Fig. 4 .Fig. 4 .

 44 Fig. 4.26 Filtered residual errors along ω = ±c 0 k for coherent wave: (a) R sd T * * W[h t h x ] , (b) R sd T S * * W[h t h x ]

  with eq. (4.48) for coherent waves (fig. 4.19) and incoherent waves (fig. 4.22). Numerical results are illustrated in fig. 4.26 for coherent waves and in fig. 4.27 for incoherent waves. Compared to the residual errors without filtering in fig. 4.20 and fig. 4.24, errors concentrated at low frequencies are largely reduced. It is observed that the errors are mainly concentrated around the main frequency band of wave fields with the effect of the filtering of window functions. It means that we get a more reasonable distribution of errors in terms of frequency.

  2.3. The numerical results are illustrated in fig. 4.28 and fig. 4.30, fig. 4.29 and fig. 4.31. It can be observed that for coherent waves (figs. 4.28 and 4.29

Fig. 4 .Fig. 4 .Fig. 4 .Fig. 4 .

 4444 Fig. 4.28 Results for coherent waves: (a) R sd T , (b) R sd T S

Fig. B. 1

 1 Fig. B.1 Example of the DFT of a periodic signal

125

 125

  Fig. B.3 uu * for point A with continuous edges: (a) one period, (b) three periods

Fig. B. 4

 4 Fig. B.4 Example of the DFT of the periodic signal presented in fig. B.1: sampled signal with a hann window function and its DFT

Re

  -i(s-τ)ω ′ dω ′ R e -i(z-y)k ′ dk ′ e i(sω+yk) dt ′ dx ′ dτdydsdz (D.1)By the formula of the Fourier transform of function f (x) = e iτx :f (s) = R e iτx e -isx dx = R e -i(s-τ)x dx = 2πδ (sτ)

Fig. E. 3 Fig. E. 4

 34 Fig. E.3 Hann window function in space and its frequency spectrum: (a) h(x), (b) | f f t(h)|

  Defining the displacement field and the Wigner transform in new coordinates by:

(F. 5 )

 5 we obtain the following relation:ũ(χ ⊥ , χ ∥ ) = u(t, x); W ε [u ε ](t, x; ω, k) = 1 c 0 Wε [ ũε ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) (F.6)

F. 1 (

 1 with the coherent forward waves, fig.4.17a) and the calculation region here is(ζ ⊥ , ζ ∥ ) = ( √ 2ζ x , √2ζx ). The Wigner transforms within both global and local coordinate systerm are compared fig. F.2. It is found that these two results are very close to each other along the characteristic line (normalized and the coefficient 1 c 0 considered). Figure F.2d shows that the numerical results calculated using the local coordinate system have less fluctuations (by spectral leakage).

Fig. F. 1 1 √ 2 Fig. F. 2 1 √ 2 ( 1 √ 2 (

 11221212 Fig. F.1 The rotational coordinates in the time-space domain (with u(t, x) of a forward wave) and in phase space (with W(t 0 , x 0 ; ω, k)):(χ ⊥ , χ ∥ ) = ( 1 √ 2 (c 0 tx), 1 √ 2 (c 0 t + x)), (k ⊥ , k ∥ ) = ( 1 √ 2 ( ω c 0k), 1 √ 2 ( ω c 0 + k))

  ∥ ) W+ [ ũε ] -W-[ ũε ]

Theorem F. 2 .

 2 The filtered residual errors in the local basis associated to the wave propagation characteristics are defined with:   Rsd T [ ũε ] * * W[ hχ ⊥ hχ ∥ ] = R+ [ ũε ] * * W[ hχ ⊥ hχ ∥ ] + R-[ ũε ] * * W[ hχ ⊥ hχ ∥ ] Rsd T S [ ũε ] * * W[ hχ ⊥ hχ ∥ ] = R+ [ ũε ] * * W[ hχ ⊥ hχ ∥ ] -R-[ ũε ] * * W[ hχ ⊥ hχ ∥ ] (F.

Fig. F. 3 Fig. F. 4 Fig. F. 5 Fig. F. 6

 3456 Fig. F.3 Residual errors along the characteristic line for coherent forward wave: (a) Rsd T , (b) Rsd T S

  High frequency setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.4 Kinetic model for high frequency wave propagation in the weak coupling limit 1.3 Error estimation for numerical solutions of wave equation . . . . . . . . . . . . . . . Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical exact errors for 1D elastic wave propagation . . . . . . . . . . . . 2.3.4 Numerical estimated errors for 1D elastic wave propagation . . . . . . . . . 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Energy-based a posteriori error estimation for high frequency wave equation in homogeneous media 3.1 Wigner transform applied to the elastic wave equation . . . . . . . . . . . . . . . . . 3.1.1 Temporal Wigner transform in signal processing . . . . . . . . . . . . . . . 3.1.2 Discussion and analysis of the interference property of Wigner transform . . 3.1.3 Scaled spatio-temporal Wigner transform in the high frequency limit . . . . . 3.1.4 Discrete Wigner transform and analysis of parameters . . . . . . . . . . . . 3.2 Transport energy-based residual errors of numerical solutions of wave equation . . . Transport equation in homogeneous media . . . . . . . . . . . . . . . . . . 3.2.2 Theoretical residual errors based on transport equation . . . . . . . . . . . . 3.3 Numerical application for wave propagation in a 1D homogeneous medium . . . . . 3.3.1 Residual errors in phase space for different refinements of finite element solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Residual errors in time-space for different refinements of finite element solutions 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Energy-based a posteriori error estimation for high frequency wave equation in heterogeneous media 4.1 Radiative transfer equation in heterogeneous media in the weak coupling limit . . . . 4.1.1 Statistical description of weak randomly-fluctuating media . . . . . . . . . .

	Table of contents List of notations and definitions 1 Introduction 1.1 Context . Table of contents 3.2.1 Introduction

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Modelling of elastic wave propagation in heterogeneous media . . . . . . . . . . . . 1.2.1 Wave propagation in homogeneous media . . . . . . . . . . . . . . . . . . . 1.2.2 Phenomena and scattering regimes for wave propagation in heterogeneous media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 1.3.1 Numerical methods for wave equation . . . . . . . . . . . . . . . . . . . . . 1.3.2 Error estimation methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 2 Explicit a posteriori error estimation for elastodynamic equation 2.1 Reconstructions in time and in space of finite element wave solutions . . . . . . . . . 2.2 A posteriori error estimate for wave equation . . . . . . . . . . . . . . . . . . . . . 2.3 Numerical application in 1D homogeneous and heterogeneous media . . . . . . . . . 2.3.1 Definition of numerical model and parameters . . . . . . . . . . . . . . . . . 2.3.2 Analytical solutions of wave equation . . . . . . . . . . . . . . . . . . . . . 2.3.3

  The second method proposes

			Introduction
	wave equation u ε	(chapter 2) classical error estimation	u hΔt FE wave solutions
	transport/radiative transfer equation W[u ε ]	residual error estimation (chapter 3,4)	Wigner transform of numerical waves W ε [u hΔt ]
	Fig. 1.6 Principle of the two a posteriori error estimations considered in this thesis

  C 1 in time and C 2 in each time interval I n

	(2.3b)
	Also, it is known that if the values of a time function and its derivatives at two extremities of each
	time interval I n are given as (w n , w n+1 , ∂ t w n , ∂ t w n+1 ), we can reconstruct a function C 1 in time, cubic
	in each time interval I n by the following equation (fig. 2.2):

Table 2

 2

		.2 Numerical parameters of media and load
	h (m)	0.55 0.28 0.14 0.07 0.035 0.017 0.008
	λ ( f c )/h	1	2	4	8	16	32	65

Table 2 .

 2 3 Mesh sizes used in the homogeneous case

	h (m)	0.066 0.035 0.017 0.008	0.004	0.002

(λ 1 ( f c )/h, λ 2 ( f c )/h) (

8,1) (16,2) (32,4) (65,8) (139,16) (278,32) Table 2.4 Mesh sizes used in both heterogeneous cases

Table 4 .

 4 2 Numerical parameters of media, load and waves for the numerical examples

	mm)

  Residual errors based on 1D radiative transfer equation in terms of W s and W d are defined as:

	4.2 Energy-based residual errors of numerical solutions of wave equation	91
	Definition 4.3.
	.40)

  ′ )+iz(k-k ′ ) dsdz dt ′ dx ′ dω ′ dk ′

	h t (t -t ′ -)e is(ω-ω = εs 2 )h x (x -x ′ -εz 2 )h * t (t -t ′ + εs 2 )h * x (x -x ′ + εz 2 1 (2π) 4 R 6 u ε (t ′ -ετ 2 , x ′ -εy 2 )u * ε (t ′ + ετ 2 , x ′ + εy 2 )h t (t -t ′ -εs 2 )h x (x -x ′ -εz 2 )
	h * t (t -t ′ +	εs 2	)h * x (x -x ′ +	εz 2	)
					1: sampled signal with a hann
	window function and its DFT

  Convolution properties of Wigner transformW ε [u ε ](t, x; ω, k) * * W ε [h t h x ](t, x; ω, k) (xx ′ -εy 2 )e i(τω+yk) dt ′ dx ′ dτdy (D.3) Changing the variables (t ′ , τ) to (α, β ) and (x ′ , y) to (m, n) with: tt ′ + ετ 2 = εα tt ′ -ετ 2 = εβ xx ′ + εy 2 = εm xx ′ -εy 2 = εn dt ′ dτ = εdαdβ , dx ′ dy = εdmdnUsing the same steps of proof above, we can also derive that:W ε [u ε ](-k) * * W ε [h t h x ] = ε 2π R 4 u ε (tεα, xεm) h * t (εα) h x (εm) e i(αω-mk) dαdm 2 = |STFT[u ε , h t h * x ](t, x; ω, -k)| 2

	=	1 (2π) 4	R 6	u ε (t ′ -	ετ 2	, x ′ -	εy 2	)u * ε (t ′ +	ετ 2	, x ′ +	εy 2	)h t (t -t ′ +	εs 2	)h x (x -x ′ +	εz 2	)
			h * t (t -t ′ +		εs 2	)h * x (x -x ′ +	εz 2	) 2πδ (s -τ) 2πδ (z -y) e i(sω+yz) dt ′ dx ′ dτdydsdz
	=	1 (2π) 2	R 4	u ε (t ′ -	ετ 2	, x ′ -	εy 2	)h * t (t -t ′ +	ετ 2	)h * x (x -x ′ +	εy 2	)
			u * ε (t ′ +	ετ 2	, x ′ +	εy 2	)h t (t -t ′ -	ετ 2	)h x (D.4)
	eq. (D.3) becomes:
			W = u = ε 2 (2π) 2 R 4 ε 2 (2π) 2 R 2 u ε (t -εα, x -εm) h * t (εα)h * x (εm) e i(αω+mk) dαdm
								u
						R 2

ε [u ε ](t, x; ω, k) * * W ε [h t h x ](t, x; ω, k) ε (tεα, xεm) h * t (εα)h * x (εm) u * ε (tεβ , xεn) h t (εβ )h x (εn) e i((α-β )ω+(m-n)k) dαdβ dmdn * ε (tεβ , xεn) h t (εβ )h x (εn) e -i(β ω+nk) dβ dn = ε 2π R 2 u ε (tεα, xεm) h * t (εα)h * x (εm) e i(αω+mk) dαdm 2 = |STFT[u ε , h t h x ](t, x; ω, k)| 2 (D.5)

  [START_REF] Babuška | A posteriori error estimates for the finite element method[END_REF] whereR+ [ ũε ] * * W[ hχ ⊥ hχ ∥ ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) = 2 √ 2Re STFT[∂ χ ∥ ũε , hχ ⊥ hχ ∥ ](χ ⊥ , χ ∥ ; -√ 2k, 0) STFT * [ ũε , hχ ⊥ hχ ∥ ](χ ⊥ , χ ∥ ; -STFT[ ũε , hχ ⊥ hχ ∥ ](χ ⊥ + εl, χ ∥εl; 0, -√ 2k) STFT * [ ũε , hχ ⊥ hχ ∥ ](χ ⊥εl, χ ∥ + εl; 0, -[ ũε ] * * W[ hχ ⊥ hχ ∥ ](χ ⊥ , χ ∥ ; k ⊥ , k ∥ ) = 2 √ 2Re STFT[∂ χ ⊥ ũε (χ ⊥ , χ ∥ ), hχ ⊥ hχ ∥ ](0, STFT[ ũε , hχ ⊥ hχ ∥ ](χ ⊥ + εl, χ ∥εl;

							√
							2k, 0)
	+ -	√ 2 16 R √ 2 16 R	STFT[ ũε , hχ ⊥	h * χ ∥ ](χ ⊥ + εl, χ ∥ -εl; -√ 2k, 0) STFT *	[ ũε , hχ ⊥	h * χ ∥ ](χ ⊥ -εl, χ ∥ + εl; -√ 2k, 0) P (2) ( √ 2l)dl √ √ 2k) P (2) ( 2l)dl (F.14)
	and
	R-√	√
				2k) STFT			2k)
	+ -	√ 2 16 R √ 2 16 R	STFT[ ũε , hχ ⊥	h * χ √ 2k) P (2) ( √ 2l)dl √ 2k, 0) STFT √ √ 2k, 0) P (2) ( 2l)dl (F.15)
	where STFT stands for STFT of ũε in new coordinates.

* [ ũε , hχ ⊥ hχ ∥ ](χ ⊥ , χ ∥ ; 0, ∥ ](χ ⊥ + εl, χ ∥εl; 0, √ 2k) STFT * [ ũε , hχ ⊥ h * χ ∥ ](χ ⊥εl, χ ∥ + εl; 0, * [ ũε , hχ ⊥ hχ ∥ ](χ ⊥εl, χ ∥ + εl;

This result is used for deducing a guaranteed upper bound of x here. In fact, x

≤ ax + b 2 ⇒ (x -a 2 ) 2 ≤ a 2 4 + b 2 ⇒ (x -a 2 ) 2 ≤ ( a 2 + b) 2 ⇒ x ⩽ a + b for a, b positive

ε in the asymptotic expansion of the Wigner measure, resulting from the magnitude of fast fluctuations of random media defined in eq. (4.1). Besides, the transport equation can be seen as a particular case of the radiative transfer equation by removing the fast fluctuations of media.First of all, we can also write the rescaled wave equation in the high frequency limit as we did in eq. (3.34):(ε∂ t ) 2 I d -Γ Γ Γ(x; ε∇ x ) • u ε (t, x) = 0 (4.6)In the case of heterogeneous media, Γ Γ Γ depends on x. However, we still present the case of isotropy for simplicity. According to the randomly fluctuating elasticity tensor defined by eq. (4.1), the Christoffel tensor Γ Γ Γ is expanded as:

Acknowledgements

, Anne-sophie Bonnet-Ben Dhia, Xavier sister for supporting me spiritually throughout finishing this thesis.

List of notations and definitions

The next list describes several notations and definitions for some spaces, products and norms that are used within the thesis.

Energy-based error estimation for HF wave in heterogeneous media

We give here the definitions of the double convolution of Wigner transform and the short-time Fourier transform: Definition 4.4. The double convolution of two spatio-temporal Wigner transforms, denoted by " * * ", is defined with respect to four variables (t, x; ω, k) by

Definition 4.5. The short-time Fourier transform is defined as [START_REF] Auger | Time-frequency toolbox[END_REF]:

The Wigner transform has the following property:

Property 8. The relation between the double convolution of Wigner transforms and the STFT is the following: ∀h(t, x),

Proof. See Appendix D.

It is obvious that we have W

The computation cost of the STFT is in the same order as that of the Wigner transform. Thus the calculation of the right-hand side of eq. (4.47) is more efficient than that of the convolution of two Wigner transforms.

If we choose w as a product of two window functions in time and in space, denoted by h(t, x) = h t (t)h x (x), the property 8 illustrates that a filtering in time-space for wave fields will filter also the results of its Wigner transform by the convolution. It means that this filtering results in smoother values in phase space. Besides, a discussion on bilinear properties of STFT can be found in Appendix E. Compared to Wigner transform, STFT suppresses largely cross terms at zero and low frequency and less frequency leakage is observed.

Based on the analysis above and by recalling that R sd T = R + + R -and R sd T S = R + -R -, the following results are derived: Theorem 4.6. Filtered residual errors based on 1D radiative transfer equation are defined as:

Appendix A

Definition of spaces and useful inequalities

In the following Ω denotes an open set of R d equipped with the Lebesgue measure x. We introduce firstly some usual definitions:

We define:

If Ω is a bound domain ⊂ R d with boundary ∂ Ω, we define:

where the space D(Ω) is the space of functions from C ∞ ( Ω) with a compact support in Ω [START_REF] Adams | Sobolev spaces[END_REF]. The function w is called the weak divergence of v and we note that

Definition of spaces and useful inequalities

Definition A.4 (H(div, Ω)). H(div, Ω) is the space of all the functions which admit the weak divergence.

Some useful inequalities are also given that are used in our work:

Theorem A.5 (Poincaré's inequality). Suppose that S an element of partition T h and h s its size, ∀u ∈ H 1 (S), there exists a constant C p such that [START_REF] Bebendorf | A note on the poincaré inequality for convex domains[END_REF][START_REF] Payne | An optimal poincaré inequality for convex domains[END_REF]:

where u S is the mean value of u in S and the norms are defined by ∥u∥ (S) = ∥u∥ L 2 (S) and ∥u∥ 1(S) = ∥∇ x • u∥ (S) . u s is introduced here since the equilibrium of mean value in each element will be used in the derivation of our error upper bound.

Theorem A.6 (Korn's inequality). Assume that Ω is a bound domain and

, u| Γ 0 = 0 , there exists a constant C k > 0 such that:

where

Theorem A.7 (Cauchy-Schwarz inequality). For all functions u 1 and u 2 in a vector space F with an inner product (, ) F , it is true that:

Appendix B

Spectral leakage in FFT

When we calculate Wigner transform in fig. 3.8, some noticeable fluctuations are found around lines that are parallel to the characteristic directions defined by ω = ±c 0 k. They are in fact the results of a well-known phenomenon called "spectral leakage", which happens in the discrete Wigner transform calculation.

The spectral leakage is a problem that arises in the digital processing of signals. In fact, for DFT, a signal is truncated firstly in order to obtain a finite set of samples. DFT implicitly assumes that the sampled signal essentially repeats itself after the measured period and hence the signal is continuous (conceptually, juxtapose the measured signal repetitively). This leads to glitches or discontinuities in the boundary of sampling period. Of course, leakage will not occur if a FFT is precisely computed on a periodic signal sampled of an integer number of cycles. However, if the measurement region is purposefully made to be a non-integral multiple of the actual signal rate, these sharp discontinuities will always spread out in the frequency domain leading to spectral leakage. In fig. B.1, when the original periodic signal with one frequency f 0 is sampled in a chosen repeating period with discontinuities in the boundaries, its unique component f 0 in frequency domain has a leakage to other frequencies.

For our wave propagation problems, the Discrete Wigner transform of u(t, x) can be seen as a 2D FT of the auto-correlation functions of displacement uu * (τ, y) in a finite area

). The yellow zone can be seen as a time-space period of FFT in 2D.

Appendix C

Pseudo-differential operator

In order to simplify the derivation of the transport equation and the radiative transfer equation with Wigner transform, pseudo-differential calculus are taken into account [START_REF] Baydoun | Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media[END_REF][START_REF] Lions | Sur les mesures de wigner[END_REF].

Definition C.1 (Pseudo-differential operators). Let ϕ(t, εD t ) and ψ(x, εD x ) be two matrix-valued pseudo-differential operators defined by:

We assume that ϕ(t, iεω) and ψ(x, iεk) are smooth functions and use the same mathematical symbols for the operators ϕ(t, εD t ), ψ(x, εD x ) and their symbols ϕ(t, iεω), ψ(x, iεk).

This definition can be seen as a more general expression of time and space differential operators. In effect, by the derivative property of the Fourier transform it is known for example that:

We have the following properties for the pseudo-differential operator [START_REF] Bal | Kinetics of scalar wave fields in random media[END_REF]:

For the heterogeneous operator ψ(x, εD x ), we have:

2) where the Poisson bracket is defined by: {ψ

, and we assume that the differential operator D x within the observable ψ acts on

2 ) should be interpreted as the component-wise inverse Fourier transform (with respect to k ′ ) of the matrix

Pseudo-differential operator If ψ(x, εD x ) = ψ(εD x ) is independent of x, the eq. (C.2) becomes:

In the same way, we can derive for ϕ(εD t ):

Appendix D

Convolution properties of Wigner transform

In this chapter, we introduce the proof of the property 8 applied in the section 4.3.

Proof. Using the definition of convolution, we have:

Properties and examples of STFT

The STFT, or the spectrogram (the magnitude squared of the STFT), is a simple extension of the FT. It can be understood as the FT is repeatedly evaluated for a running windowed version of the time-space domain signal. Cross terms exist also for STFT. Consider a general case of the sum of two signals u 1 + u 2 , its STFT has the following form:

where |STFT[u i , h tx ]| and φ u i (t, x; ω, k), i = 1, 2 are respectively the magnitude and the phase of the STFT of signal u i . In eq. (E.1), the first two terms are the auto-terms and the last term is the cross term. It can be observed that if two STFTs have no overlapping (t, x)-(ω, k) support, the cross term is zero. Equation (E.1) can also be generalized to sum of N signals.

A summary of the principal properties of cross terms of STFT are given here. More detailed analysis can be found in [START_REF] Kadambe | A comparison of the existence of'cross terms' in the wigner distribution and the squared magnitude of the wavelet transform and the short-time fourier transform[END_REF][START_REF] Lu | Deconvolutive short-time fourier transform spectrogram[END_REF][START_REF] Sandsten | Time-frequency analysis of time-varying signals and non-stationary processes[END_REF].

1. The STFT cross terms occurs at the intersection of two overlapping transforms and no cross term occurs between two frequency components like Wigner transform in eq. (3.13). Depending upon the amount of overlap of each pair of transforms, the energy distribution of the STFT of an N component signal can have a minimum of zeros and a maximum of N(N-1) 2 cross terms, unlike the Wigner transform, which always has N(N-1) 2 cross terms in the midway of two Wigner transform auto-terms. However, it also reduces the time-frequency resolution simultaneously.

2. The STFT cross terms can have a maximum magnitude equal to twice the product of the magnitude of the two spectrograms.

3. The STFT cross terms are modulated by a cosine whose argument is a function of the difference of phases of STFT[u i , h tx ].

Properties and examples of STFT

Here we study the same example, namely the sum of two complex-valued and separate gaussian signals F(t), as we presented for the Wigner transform in the section 3.1.2 to view the cross term properties of STFT or spectrogram (hann window is used here with length = sample length/6). Properties and examples of STFT is due to superposition of cross terms). However, it is known that it cannot suppress all the cross terms between positive frequencies because a continuous frequency range is used here. Finally, less frequency leakage is observed for spectrogram due to the use of window functions. Wigner transform is calculated directly in a limited calculation region and we can treat it as application of a rectangular window, whose frequency spreading is very wide, which may mask important spectrum details at even lower levels. That's why we find small fluctuations everywhere in fig. E.4b. Generally, non-rectangular window functions actually increase the total leakage, but they can also redistribute it to places where it does the least harm. To different degrees they reduce the level of the spreading by increasing the high-level leakage in the near vicinity of the original component [START_REF] Harris | On the use of windows for harmonic analysis with the discrete fourier transform[END_REF]. So in fig. E.4c, the use of hann window function increases leakage near main lobe and also allows to reduce frequency leakage at a remote area of the main lobe.

Appendix F

Energy-based residual errors calculated within the local basis defined by the wave propagation characteristics

In this section, we consider the residual errors in 1D heterogeneous media defined in chapter 4 using the local basis defined by the wave propagation characteristics. A transformation of coordinates is firstly performed for the Wigner transform and the radiative transfer equation in 1D case. Then the residual errors are re-calculated respectively for coherent waves and incoherent waves.

We propose a change of variables for (t, x) → (χ ⊥ , χ ∥ ) (fig. F.1):

and for (τ, y) → (φ , ψ):

where the subscripts ∥ and ⊥ represent respectively the directions along the forward and backward wave propagation characteristics. Recall that τ and y are the shift variables related to t and x in the definition of Wigner transform.

Introducing eq. (F.1) and eq. (F.2) into the definition 3.3 of Wigner transform in 1D:

Appendix G

Résumé substantiel

La simulation numérique de la propagation d'ondes en milieu hétérogène est un problème difficile car il nécessite de suivre les différents fronts d'ondes diffractés sur les hétérogénéités. C'est particulièrement vrai en régime mésoscopique à haute fréquence, car l'interaction des ondes avec le milieu est plus forte et les distances de propagation (relativement à la longueur d'onde) sont plus importantes. Ces travaux de thèse portent sur le développement d'outils d'estimation d'erreurs pour les solutions numériques d'équations d'ondes par éléments finis. Ces outils doivent permettre de contrôler la qualité des simulations et de piloter l'optimisation du coût numérique.

Dans une première partie, un estimateur d'erreur explicite a posteriori est proposé. Il permet d'évaluer une erreur en norme énergétique en espace et en norme L ∞ en temps entre la solution exacte de l'équation élastodynamique et une solution approchée reconstruite. L'idée principale du développement est basée sur une erreur en résidu et l'utilisation de reconstructions des champs. Pour tenir compte des régularités nécessaires à l'obtention d'une borne supérieure de l'erreur, plusieurs reconstructions en temps pour le déplacement et la vitesse sont proposées, ainsi que des reconstructions en espace-temps pour la contrainte. La performance de l'estimateur a été comparée par rapport à des erreurs numériques exactes pour des cas homogène et hétérogène, et pour des solutions éléments finis obtenues avec différents tailles de maillage et pas de temps. L'erreur estimée diminue bien avec le raffinement en temps et en espace, mais l'efficacité (rapport entre l'erreur estimée et l'erreur réelle) augmente également avec le raffinement. La convergence asymptotique de l'estimateur obtenu doit donc être améliorée. L'analyse des résultats montre que les composantes des erreurs liées aux résidus de l'équation d'équilibre dérivée en temps sont dominantes dans la contribution des erreurs totales estimées. show that it gives rise to a fully computable upper bound. However, its effectivity index and its asymptotic accuracy remain to be improved.

The second error estimator is derived for high frequency wave propagation problem in heterogeneous media in the weak coupling regime. It is a new residual-type error based on the radiative transfer equation, which is derived by a multi-scale asymptotic expansion of the wave equation in terms of the spatio-temporal Wigner transforms of wave fields. The residual errors are in terms of angularly resolved energy quantities of numerical solutions of waves by finite element method. Numerical calculations of the defined errors in 1D homogeneous and heterogeneous media allow validating the proposed error estimation approach. The application field of this work is the numerical modelling of the seismic wave propagation in geophysical media or the ultrasonic wave propagation in polycrystalline materials.