L Campoli

P Quemar
email: quemar@math.univ-paris13.fr

A Bonfiglioli
email: aldo.bonfiglioli@unibas.it

M Ricchiuto
email: mario.ricchiuto@inria.fr

Shock-fitting and predictor-corrector explicit ALE Residual Distribution

Introduction and main results

This work is part of a long term effort to develop a modern unstructured shockfitting algorithm. The authors have so far shown the capabilities of their approach for steady flows in two [START_REF] Paciorri | A shock-fitting technique for 2D unstructured grids[END_REF][START_REF] Paciorri | Shock interaction computations on unstructured, twodimensional grids using a shock-fitting technique[END_REF] and three [START_REF] Bonfiglioli | An unstructured, three-dimensional, shock-fitting solver for hypersonic flows[END_REF] spatial dimensions, as well as some potential for time dependent moving shocks [START_REF] Bonfiglioli | Unsteady shock-fitting for unstructured grids[END_REF]. In all the mentioned works fitted shocks are treated as interior boundaries of zero thickness that are free to move throughout a triangular/tetrahedral mesh. This underlying unstructured triangulation/tetrahedrization covers the entire computational domain, and locally adapts to exactly fit the shock motion. The latter is ruled by the Rankine-Hugoniot jump relations which are in turn deduced from an approximate flow solution obtained by discretizing the compressible flow equations on either side of the discontinuity. To this end, in this work we will focus on the numerical solution of the perfect gas Euler equations. This contribution is a first step to evaluate the unsteady fitting method when schemes with different properties are used. We evaluate the results obtained with explicit residual distribution schemes in Arbitrary Lagrangian Eulerian (ALE) form, developed in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Arpaia | An ALE formulation for explicit rungekutta residual distribution[END_REF].

The paper investigates quantitatively and qualitatively the influence on the fitted solutions of the accuracy order (first or second here) of the method, as well as on the stabilization mechanism used. The main result is that accuracy seems to be the most influent factor.

2 Numerical methods

Shock-fitting algorithm

The unstructured shock-fitting algorithm consists of two key ingredients: 1) a local re-meshing technique that constructs a time-dependent mesh in which the fitted discontinuities are internal boundaries of zero thickness and 2) an algorithm for solving the Rankine-Hugoniot jump relations that provide the Lagrangian velocity of the discontinuity and an updated set of dependent variables within the downstream side of the fitted shock. More precisely, in two space dimensions, the fitted shock fronts are made of polygonal curves, i.e., a connected series of line segments (which we call the shock edges) that join the shock points. These shocks are free to move throughout a background triangular mesh that covers the entire computational domain (see Fig. 1a). Starting from a background grid, at a given time level n local mesh modifications are performed by re-generating a constrained Delaunay triangulation in the neighbourhood of the shock. This allows to ensure that the edges that make up the shock front are also part of the mesh that covers the entire computational domain (see Fig. 1b). The shock speed σ sh is computed from the Rankine-Hugoniot jump relations within each of the shock points. Using nodal values of shock speed and shock normal n, the fitted shock front is moved in a Lagrangian manner, as shown in Fig. 1c. A second-order-accurate temporal integration of the shock trajectory is obtained using a predictor-corrector scheme [START_REF] Bonfiglioli | Unsteady shock-fitting for unstructured grids[END_REF].

In the smooth regions on either side of the shock, a vertex-centred ALE discretization, described in Section 2.2, allows to obtain an approximation of the flow variables on the mesh. The ALE formulation allows to embed naturally in the schemes the presence of the shock front as a moving boundary. This is schematically shown in Fig. 1c where dashed lines are used to show the edges of the triangular cells when the shock front is at time level n and solid lines to show the same sides when the shock has reached time level n + 1.

The interested reader can refer to [START_REF] Bonfiglioli | Unsteady shock-fitting for unstructured grids[END_REF] for additional details on the unsteady fitting procedure.

Predictor-corrector ALE residual distribution

We consider discretization methods for the Euler equations recast in the following ALE compact form:

∂ t (Jw) + J∇ • (f(w) -σw) = 0 (1)
with w the array of conserved quantities (mass, momentum and energy), f the conservative fluxes, J the determinant of the Jacobian of the transformation between the reference and actual frame, and σ the local deformation velocity. To discretize (1), we consider the two-step explicit Residual Distribution (RD) method developed in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Arpaia | An ale formulation for explicit runge-kutta residual distribution[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF],which can be written in the compact discrete form

|C n+1 i |w n+1 i =|C n+1 i |w * i -∆t ∑ K∈D i Φ K i (w n h , w * h ) (2)
where the fluctuations Φ K i define a splitting of the average element residual Φ K defined by

∑ j∈K Φ K j = Φ K = 1 ∆t K n+1 w * h - K n w n h + 1 2 Φ K (w n h ) + 1 2 Φ K (w * h ) (3)
with

∑ j∈K Φ K j = Φ K = ∂ K n+1/2 f(w h ) -σ h w h • n ds (4)
The vector σ h represents the mesh velocity with nodal values σ i = (x n+1 i x n i )/∆t . The values w * i are obtained from a first order predictor which is computed as

|C n+1 i |w * i = |C n+1 i |w n i -∆t ∑ K∈D i ΦK i (w n h ) (5)
where now the fluctuations ΦK i are a splitting of the following geometrically nonconservative average steady cell residual [START_REF] Arpaia | An ale formulation for explicit runge-kutta residual distribution[END_REF] ∑ j∈K

ΦK j = ΦK = ∂ K n+1/2 f(w h ) • n ds - K n+1/2 σ h • ∇w h dx (6)
Almost all the geometrical quantities in the scheme are evaluated on the half-time n + 1/2 averaged configuration T n+1/2 h which ensures

|K n+1 | -|K n | = ∆t ∂ K n+1/2 σ h • n ds = 0 , (7)
thus the satisfaction of a Discrete Geometric Conservation Law (DGCL).

The properties of the method introduced above are determined by the definition of the split residuals Φ K i , Φ K i and ΦK i , which, following [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF][START_REF] Arpaia | An ale formulation for explicit runge-kutta residual distribution[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF] have the form :

Φ K i = ∑ j∈K m K n+1 i j w * j -m K n i j w n j ∆t + 1 2 Φ K i (w n h ) + 1 2 Φ K i (w * h ) with Φ K i = β i Φ K (w n h ) (8)
and similarly for ΦK i . The β j are distribution matrices, uniformly bounded w.r.t the cell residuals (6), (3), while the m K i j s are mass matrix entries consistent with the definition of the spatial distribution [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF]. Independently of the specific construction of these quantities, the above definitions give a scheme which is formally second order accurate in space and time, fully conservative, and verifying the DGCL. The interested reader is referred to [START_REF] Arpaia | An ale formulation for explicit runge-kutta residual distribution[END_REF][START_REF] Arpaia | r-adaptation for Shallow Water flows: conservation, well balancedness, efficiency[END_REF][START_REF] Ricchiuto | An explicit residual based approach for shallow water flows[END_REF][START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF] for more details. In practice, we have used two main classes of methods: multidimensional upwind (MU) schemes, and centered schemes with some form of stabilization. In particular, the computations shown in the paper will compare results using the first order linear and monotone N scheme and the second order linear LDA scheme. These are MU distribution methods. A nonlinear method obtained by blending these two, and referred to as the LDAN scheme, is also tested. Two non-upwind methods are also considered, in particular, the explicit predictor-corrector formulation of the second order linear Streamline-Upwind (SU) method proposed in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF] (see also [START_REF] Hughes | Stabilized methods for compressible flows[END_REF] and references therein), and the nonlinear blended central (Bc) discretization obtained when blending the SU method with a limited Lax-Friedrich's distribution. The interested reader may consult the above mentioned references for more details.

3 Numerical results

Shock-expansion interaction

We start by investigating the behaviour of the schemes on the computation of an expansion fan interacting with a moving shock. The problem involves a rectangular cavity initially filled with an inviscid perfect gas at rest, as shown in the sketch of Fig. 2a. The left end of the domain acts as a piston which is impulsively set into motion, generating a planar shock wave moving towards the closed (right) end of the domain. Given the density ratio across the shock wave, all other kinematic and thermodynamic quantities follow from the Rankine-Hugoniot (R-H) jump relations and are summarized in Tab. 1. The upstream and downstream dimensionless quantities have been labelled with the subscripts 1 and 2 in Tab. 1. The Mach number in a reference frame attached to the shock is denoted by M r , and its upstream value coincides with that of the dimensionless shock speed W s , since a 1 has been chosen as the reference velocity.

After the impulsive startup, the piston starts slowing down at t = t 0 + t s ; its velocity and position are prescribed analytically by:

u p (t) = u 2 t 0 < t ≤ 0 u 2 exp -t t ≥ 0 , x p (t) = u 2 t t 0 < t ≤ 0 u 2 (1 -exp -t ) t ≥ 0 , (9)
where t 0 = -t s and:

t s = x 0 sh W s -u 2 0.3755 ρ u v p M r a Upstream (1)
1.4 0.0 0.0 1.0 1.5336 1.0 Downstream (2) 2.6872 0.7346 0.0 2.5772 0.6895 1.1588 Table 1: Shock-expansion interaction: initial conditions.

The flow evolution is sketched in the xt plane of Fig. 2c. We start the simulations at t = 0 with the shock located x 0 sh = 0.3 unit lenghts ahead of the piston and we initialize the solution using the initial conditions listed in Tab. 1.

The tail of the expansion wave (marked by arrows in Fig. 2c) is the straight, characteristic line of slope u 2 + a 2 that is shed from the origin of the space-time plane and divides the expansion wave from the uniform, shock-downstream region, labelled 2 in Fig. 2. Up to time:

t = x 0 sh (u 2 + a 2 ) -W = x 0 sh (1 -M r,2 ) a 2 0.8339
when the tail of the expansion wave meets the shock trajectory, the expansion region is a simple wave, so that an analytical solution can be computed as described, for instance, in [20, §4.1.5]. At later times, the expansion wave is reflected off the moving shock as an expansion wave of the opposite family, the shock starts reducing its speed and the region of the xt plane that is between the moving piston and the moving shock ceases to be a simple wave.

Simulations have been advanced until the final time t 2 = 2.8 when the shock is located at x t 2 sh 4.43 and has almost reached the closed (right) end of the cylinder. We present results obtained on a background grid of 639 nodes, 12448 triangles, with a shock front discretized by 41 points (cfr. right picture of Fig. 2a). During the time evolution, the streamwise size of the computational domain is reduced because of the piston's motion. All nodes of the background mesh move as the time elapses, but the topology of the mesh does not change. At any time instant, the location of the nodes of the background mesh is prescribed using the following equation:

x i (t) = x i (0) + x p (t) 1 - x i (0) L
where L = 5 is the streamwise length of the computational domain at t = 0. Figure 3 provides the pressure field at three subsequent time instants for all the five different spatial discretization schemes and the two shock-modelling schemes. Apart from the excessive shock-thickness of the shock-capturing solutions, it is clear that shock-fitting provides a much cleaner representation of the shock-downstream region. Observe, in particular, that the crossflow non-uniformities are considerably less pronounced in the shock-fitting calculations.

Figure 4 shows the streamwise pressure profiles at two subsequent time instants for all the five spatial discretization schemes. Note that the figures report all the data points in the mesh. The superiority of the shock-fitting calculations is evident, particularly in the reduction of cross-shock spurious variations, despite the fact that both the fitted, and captured solutions are represented on unstructured triangulations with almost identical topologies. We do observe an impact of the numerical discretization chosen. The upwind or centered nature of the method does not seem to affect the results as much as the resolution of the scheme. The N scheme clearly provides the "wrong boundary conditions" for the shock, which seems slighly miscplaced w.r.t. the analytical solution, as one can see in the right frame of Fig. 4a. Shock amplitude and position are instead perfectly matched by the high order schemes, with perhaps a barely visible off-set present only in the LDAN result, shown in Fig. 4c. These differences are clearly much more important in the captured solutions: the linear, first-order-accurate N scheme smears the shock-profile over a considerable streamwise extent; as expected, the linear, secondorder-accurate, but non-monotone LDA and SU schemes give rise to overshoots around the discontinuity; the non-linear, second-order-accurate LDAN scheme is monotonicity preserving, but the otherwise similar Bc scheme also overshoots the shock profile.

Shock-vortex interaction

The interaction between a shock and a vortex has been frequently reported in the literature as a tool for understanding the mechanisms of noise generation due to the interaction between a shock-wave and a turbulent flow [START_REF] Pao | A numerical study of two-dimensional shock vortex interaction[END_REF][START_REF] Meadows | Computational study on the interaction between a vortex and a shock wave[END_REF][START_REF] Inoue | Sound generation by shock-vortex interactions[END_REF][START_REF] Grasso | Shock-wave-vortex interactions: Shock and vortex deformations, and sound production[END_REF]. This problem has also been frequently used as a code verification benchmark [23, [START_REF] Abgrall | Construction of second order accurate monotone and stable residual distribution schemes for unsteady flow problems[END_REF][START_REF] Rossiello | Second-order-accurate explicit fluctuation splitting schemes for unsteady problems[END_REF]. Shock-fitting calculations of shock-vortex interactions have been previously reported in [START_REF] Pao | A numerical study of two-dimensional shock vortex interaction[END_REF][START_REF] Salas | Shock-fitted Euler solutions to shock-vortex interactions[END_REF][START_REF] Arina | A shock-fitting spectral method for the study of shock-vortex interactions[END_REF].

A uniform, supersonic stream, characterised by a shock-upstream Mach number M s = |u ∞ |/a ∞ , carries a vortex, from the left to the right of Fig. 5, towards a stationary normal shock. The computational domain, which is sketched in Fig. 5 along with the boundary conditions applied along its boundaries, is the rectangle [0, 2L] × [0, L]. At the initial time, t = 0, the vortex is centred in (x v /L, y v /L) = (0.5, 0.5) and the shock is located 0.2 L downstream of the vortex centre. The vortical structure we have used, which is an exact solution of the steady Euler equations, consists in a clockwise vortex characterised by a purely tangential velocity component. More precisely, using a polar coordinate system with the origin at the centre of the vortex and moving at constant speed |u ∞ |, the perturbation velocity field ũ reads:

ũθ = -ε |u ∞ | τ e α(1-τ 2 ) (10a) ũr = 0 (10b)
In Eq. ( 10) τ = r/r c is the non-dimensional radial distance from the pole of the moving reference frame and ε, α and r c are user-defined parameters that control the shape and magnitude of the perturbation. This particular solution of the Euler equations features a divergence-free velocity field, which implies that density is constant along the streamlines and a non zero vorticity field. Since the free-stream flow is isentropic, the thermodynamic variables can be easily obtained from the linear momentum equation, and Crocco's form of the steady momentum equation implies that there must be a gradient of the perturbation total enthalpy in the radial direction, so that total enthalpy also changes across the streamlines.

A vortex Mach number, M v , can be defined as the ratio between the maximum velocity perturbation and the sound speed, a ∞ , of the shock-upstream flow: M v = max ( ũθ ) /a ∞ . For the particular vortical structure used here, the vortex Mach number becomes:

M v = ε M s e (α-1 2 ) √ 2 α (11)
Making the same choice of constants in Eqs. ( 10) and ( 11) also used in [23, 1, 18]: α = 0.204 and r c /L = 0.05, the vortex radius is about 0.35 L, the maximum velocity perturbation occurs at about 0.11 L from the vortex centre and the strength ε of the velocity perturbation follows from Eq. ( 11), once M s and M v have been set. The topological pattern that arises once the vortex impinges on the shock depends upon both the shock and vortex strengths. We shall hereafter refer to the taxonomy adopted by Grasso and Pirozzoli [START_REF] Grasso | Shock-wave-vortex interactions: Shock and vortex deformations, and sound production[END_REF], who define weak shock-vortex interactions as those that do not exhibit any shock reflection and strong shock-vortex interactions as those that feature secondary shocks; strong interactions can be further classified depending on whether a regular reflection or a Mach reflection occurs. The same authors identify in the (M s , M v ) pair the two independent variables that govern the interaction: when the vortex Mach number is sufficiently low, weak interactions always occur; however, at a given shock strength, an increase in the vortex Mach number above the M v ≈ 0.1 ÷ 0.2 threshold causes the shock to fold thus producing reflected and diffracted shocks that yield either a regular or a Mach reflection, depending upon the value of M v .

In the present work, we have set the shock and vortex Mach numbers respectively equal to M s = 2.0 and M v = 0.2, so that a weak interaction occurs. The reason for doing so is that the present shock-fitting algorithm is currently unable to automatically identify and follow the appearance of those secondary shocks that arise when a strong interaction takes place. Nonetheless, as demonstrated in [START_REF] Bonfiglioli | Unsteady shock-fitting for unstructured grids[END_REF], strong interactions can be simulated with the available algorithm, thanks to a hybrid mode of operation whereby the normal shock is fitted and the secondary shocks are captured. This hybrid modelling is however not suited for the grid-convergence analysis which will be performed hereafter, since it is impossible to separately identify the effects of the two different shock-modelling practices. Therefore, only weak interactions will be considered here and simulations will be performed using each of the five schemes, both in shock-fitting and shock-capturing mode.

When simulations are run in shock-fitting mode, the linear LDA and SU schemes, that are not monotonicity preserving, can be used without giving rise to oscillations because they are used in smooth regions of the flowfield. The use of non-linear schemes is unnecessary in this case, but non-linear scheme have also been tested in conjunction with shock-fitting to allow a comparison with the shock-capturing calculation using the same scheme.

When simulations are run in shock-capturing mode, it should be expected and it is indeed observed that the linear schemes give rise to wiggles around the captured shock wave. Nonetheless, shock-capturing calculations with linear schemes have also been performed to allow a comparison with the shock-fitting calculation using the same scheme.

A qualitative view of the solutions obtained with the two approaches is given in Fig. 6. The pictures show the total enthalpy contours in the solutions obtained with the LDA and LDAN schemes with shock capturing and fitting on the finest mesh level of table 2. Note that the fitted results are virtually indistinguishable, and only those of the LDA are reported. As for the shock-expansion interaction, the qualitative difference is already quite striking from the comparison of the contours. In particular, besides the oscillations related to the approximation of the shock, we can see clearly that the contours downstream of the discontinuity are much less smooth in the captured solutions. This is more visible in the LDAN results (Figs. 6d-6f), however a similar behaviour is also present in the LDA solutions. The fitted computations, on the other hand, show very nice and smooth contours.

In order to perform a quantitative comparison among the various combinations of spatial discretization schemes and shock-modelling practices, we have used the following global measure:

S h (t) = 1 |Ω | Ω s h (x,t) dΩ (12)
where s = pρ -γ is entropy and we have computed the corresponding discretization error and order-of-convergence. To perform a grid convergence analysis, the computational domain has been discretized using a sequence of four nested Delaunay triangulations, hereafter labelled g 1 to g 4 . The coarsest g 1 mesh has been created using the Triangle code [START_REF] Shewchuk | Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator[END_REF] by specifying a maximum area constraint, whereas the other three grid levels have been obtained by recursive subdivision of each triangle of the coarser level into four sub-triangles. Table 2 summarizes the characteristics of the various meshes that have been used both in the shock-capturing computation and as background triangulations in the shock-fitting computation. The dimensionless, characteristic spatial size h k of grid k has been set equal to the uniform spacing along the boundaries of the computational domain; observe that both h k and the dimensionless time-step length ∆t k are halved when moving from grid level k to the finer grid level k + 1.

Whenever an exact, or manufactured solution to the set of governing PDEs is available, the discretization error can be evaluated straightforwardly. However, when the exact solution is not known, the discretization error of a numerical solution and the convergence properties of a numerical scheme can still be estimated using a suite of techniques collectively known as Richardson extrapolation (RE). There exist different variations of the RE technique depending on the number of grid levels involved and on the expression of the leading error terms that is a priori assumed in the analysis; hereafter we refer to the generalized RE technque of fixed order, following the nomenclature used in [19]. The term generalized refers to the fact that the RE technique is applied to nth-order-accurate schemes and grids refined by an arbitrary factor r. Note that the use of nested triangulations fixes the grid refinement ratio between each couple of meshes in the sequence to a value r = 2. It is also assumed that the order of convergence n is known (hence, the term fixed) and is typically set equal to the design order of the scheme. The method used here relies on an error expansion of the type

ε h (x) = u h (x) -u 0 (x) = ∞ ∑ p=n g p (x) h p = g n (x) h n + HOT (13)
Neglecting the HOT , writing Eq. ( 13) for two grid levels h and rh, and solving for the unknown exact solution, one obtains:

ũ0 = u h - u rh -u h (r n -1) (14)
In Eq. ( 14) we have used a tilde to underline the fact that ũ0 is only an approximation of the corresponding quantity that appears in Eq. ( 13), because the higher-order contributions have been neglected. Since an exact solution is not known, for the present test-case, Eq. ( 14) has been used as follows: S0 = S h 4 -S h 3 -S h 4 3 to provide an approximation of the exact entropy integral of Eq. ( 12). This is a rather common practice that, according to Roache [START_REF] Roache | Quantification of uncertainty in computational fluid dynamics[END_REF], has been in use in the CFD arena since the early 1990s, even though the procedure is somewhat questionable, because the order is a-priori given in order to perform the RE.

The results obtained with this procedure are summarized in the convergence plots of Fig. 7. The plots show that all the captured solutions converge at best with a first order rate. When removing the error due to the shock-capturing, all the schemes behave as if the solution were smooth. In particular, all the schemes provide convergence rates very close to those observed in [START_REF] Ricchiuto | Explicit Runge-Kutta residual distribution schemes for time dependent problems: Second order case[END_REF]: 2 for the liner LDA and SU schemes, very close to 2 for the Bc scheme, a rate of roughly 1.5 for the LDAN scheme, and an order 1 for the N scheme.

Conclusion and outlook

We have provided a preliminary study of the behaviour of explicit ALE residual distribution when coupled with shock-fitting. The results show that the fitted solutions are quite independent on the stabilization mechanism used. Multidimensional upwind schemes, and centered schemes with some form of artificial dissipation yield very similar results. What makes most of the difference is the accuracy of the discretization. Even with fitting, first order schemes fail to provide the correct boundary conditions for the discontinuity, producing a visible error in shock position and strength. high order discretizations, both linear and nonlinear, are thus required. A comparison with captured solutions, shows two main results: fitting allows to recover second order of accuracy where capturing leads to first order; the fitted results are much less sensitive to the use of monotonicity preserving methods, as well as to mesh irregularities due to the absence of mesh dependent errors generated in the captured region. In addition anomalies related to the mesh dependent nature of the capturing error error are removed so that the methods behave as if the flow were smooth. This has been proven on two test cases.

Natural extensions of this work will be to consider more complex cases, and more interesting discretization strategies. The first aspect may involve stronger shockvortex interactions, with the appearance of secondary (lambda) shocks. Concerning the discretization, the ALE framework used here would allow quite naturally the use of a hybrid fitting/adaptive approach, with an r-adaptation strategy in the nonfitted domains, as proposed e.g. in [START_REF] Arpaia | r-adaptation for Shallow Water flows: conservation, well balancedness, efficiency[END_REF]. Lastly, a very interesting test would be to repeat the studies of this paper when the fitting is performed using non-conservative formulations in the CFD part. This may reveal very useful when considering flows also involving strong contact discontinuities, and chemistry.

Fig. 1 :

 1 Fig. 1: Shock-fitting: schematic illustration of some of the algorithmic ingredients.

  (a) Computational set-up at t = 0. (b) Nested meshes and shock position at t = 2.8.(c) Pressure iso-contours in the space-time plane.

Fig. 2 :

 2 Fig. 2: Shock-expansion interaction.

Fig. 3 :

 3 Fig.3: Shock-expansion interaction: pressure contours at times t = 0.4 (left), t = 1.2 (centre) and t = 2.0 (right): 20 equally spaced contour levels between 0.5 and 2.6, fitted solutions on top and captured solutions on the bottom of each frame.

Fig. 4 :

 4 Fig.4: Shock-expansion interaction: streamwise pressure profiles at two time instants for all the five spatial discretization schemes. Right column: zoom of the t = 0.8 peak.

Fig. 5 :

 5 Fig. 5: Shock-vortex interaction: computational domain and boundary conditions.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Shock-vortex interaction: total enthalpy contours at times t = 0.2 (left), t = 0.3 (centre) and t = 0.4 (right): 15 equally spaced contour levels obtained with the shock-capturing versions of the LDA (top), and LDAN (middle), and with shockfitting and LDA (bottom).

Table 2 :

 2 Shock-vortex interaction: grid features.

	grid N. of triangles N. of gridpoints ∆t k	h k
	g 1	6794	4181	0.00160 0.03000
	g 2	27128	14437	0.00080 0.01500
	g 3	108438	55277	0.00040 0.00750
	g 4	433778	218269	0.00020 0.00375

Acknowledgements The authors are indebted to R. Paciorri for the many constructive discussions during the development of this work.


Imen Nasr
email: nasr.imen@supcom.rnu.tn

Leïla Najjar Atallah

Sofiane Cherif
email: sofiane.cherif@supcom.rnu.tn

Benoît Geller
email: benoit.geller@ensta-paristech.fr

Hybrid Cramer Rao Bound on Time Delay Estimation in Rayleigh Fading Channels

We propose to derive the analytical expression of the Hybrid Cramer-Rao Bound HCRB for the joint time delay and channel fading estimation. The time delay is considered as a deterministic parameter and the channel is considered as a complex gaussian random variable. The HCRB expression is computed for both the data-aided (DA) and the non dataaided (NDA) modes. The obtained result allows us to assess the performance of on-line and off-line dynamical Rayleigh channel gains estimation. It also evaluates timing recovery techniques under a time varying fading channel.

I. INTRODUCTION

The standard Cramer-Rao Bound (CRB) [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF] and the Bayesian Cramer-Rao Bound (BCRB) [START_REF] Muravchik | Posterior Cramer-Rao bounds for discrete-time nonlinear filtering[END_REF] have been introduced in the literature to find a theoretical minimum mean square error (MSE) limit for respectively, deterministic and stochastic parameter estimation. Closed form expressions of the CRB have been derived in [START_REF] Bellili | Closed-form CRLBs for CFO and phase estimation from turbo-coded square-qam-modulated transmissions[END_REF] for carrier frequency and phase offset estimation for turbo-coded Square-QAM modulated signals. It was also evaluated in [START_REF]Closed-form Cramer Rao lower bound for SNR estimation from turbo-coded BPSK-, MSK-, and square-QAM-modulated signals[END_REF] and [START_REF] Bellili | Cramer-Rao Lower Bounds for Frequency and Phase NDA Estimation From Arbitrary Square QAM-Modulated Signals[END_REF] respectively for SNR and for frequency and phase NDA estimation. For the timing recovery problem, the CRB has been evaluated for DA [START_REF] Jesupret | Digital demodulator synchronization -Performance analysis[END_REF], NDA [START_REF] Masmoudi | Closed-form expressions for the exact Cramer-Rao bounds of timing recovery estimators from BPSK, MSK and square-QAM transmissions[END_REF] and with the help of soft information [START_REF] Geller | Block turbo codes: From architecture to application[END_REF], [START_REF] Diatta | Reed Solomon turbo codes for high data rate transmission[END_REF] in code-aided (CA) [START_REF] Nasr | Performance study of a near maximum likelihood code-aided timing recovery technique[END_REF], [START_REF] Nasr | CRB derivation and new code-aided timing recovery technique for QAM modulated signals[END_REF] only for a constant delay. The BCRB has been derived for the unknown random phase offset problem [START_REF] Geller | Advanced synchronization techniques for the internet of things[END_REF], [START_REF] Geller | Proc. of OCEANS IEEE Oceanic engineering, title=Equalizer for High Rate Transmission[END_REF] in [START_REF] Yang | Bayesian and hybrid Cramer-Rao bounds for the carrier recovery under dynamic phase uncertain channels[END_REF]- [START_REF] Yang | Bayesian and hybrid Cramer-Rao bounds for QAM dynamical phase estimation[END_REF] and for dynamical timing recovery [START_REF] Nasr | A soft maximum likelihood technique for time delay recovery[END_REF], [START_REF] Nasr | Near MAP Dynamical Delay Estimator and Bayesian CRB for Coded QAM Signals[END_REF].

However, it is frequent that a mixture of deterministic and stochastic parameters has to be jointly estimated. For instance, in real systems, time delay recovery is often associated to dynamic channel estimation. CRB and BCRB are then no more relevant for such estimators designed to jointly recover deterministic and random parameters. By the introduction of a general framework to derive the analytical expression of the BCRB suited for a vector of parameters estimation [START_REF] Muravchik | Posterior Cramer-Rao bounds for discrete-time nonlinear filtering[END_REF], an extension has been proposed [START_REF] Rockah | Array shape calibration using sources in unknown locations-part i: Far-field sources[END_REF] for some practical cases where the vector to be estimated is composed of a mixture of random and deterministic parts. This theoretical bound is called the Hybrid CRB (HCRB).

In [START_REF] Hijazi | Analytical Analysis of Bayesian Cramer-Rao Bound for Dynamical Rayleigh Channel Complex Gains Estimation in OFDM System[END_REF] a HCRB has been introduced for a dynamical time varying channel, modeled by a random variable, and a deterministic carrier frequency offset estimation in OFDM systems. In [START_REF] Yang | Bayesian and hybrid Cramer-Rao bounds for the carrier recovery under dynamic phase uncertain channels[END_REF], [START_REF] Bay | On the hybrid Cramer-Rao bound and its application to dynamical phase estimation[END_REF], this bound has been computed for a dynamical phase offset estimation following a brownian model with a linear unknown drift. The authors of [START_REF] Vilà-Valls | On-line Hybrid Cramér-Rao Bound for oversampled dynamical phase and frequency offset estimation[END_REF] proposed a HCRB for the joint estimation of the pair of dynamical carrier phase-Doppler shift (random part) and the constant time delay (deterministic part). In this proposal, we provide an analytical derivation of the HCRB for a joint estimation of a Rayleigh fading channel and a deterministic time delay.

This paper is organized as follows. In section II, the system model is presented. In section III, closed form expressions of the HCRB for both the DA and the NDA modes are derived. Simulation results are provided in section IV and validate our analysis. The last section concludes our work.

GLOSSARY OF PRINCIPAL NOTATIONS

• ℜ{z} and ℑ{z}: the real and imaginary parts of z • z * : the conjugate of z • X T : the transpose of the matrix X • ḟ (resp. f ): the first (resp. the second) derivative of f , for any function f • f X (x): the probability density function (pdf) of the random variable X • P (X): the probability of the random variable X • Ey(X): the mean value of the random variable X over y • |X|: the determinant of the matrix X • ∇ x : the gradient with respect to x • △ y x = ∇x∇y: the Laplace operator with repect to yand x

II. SYSTEM MODEL

Let us consider the linearly modulated transmitted signal s(t) such as:

s(t) = i a i h(t -iT ), (1)
where a i denotes the transmitted symbols which are assumed to be statistically independent and equally likely, with normalized energy, h(t) is the impulse response of the root Nyquist transmission filter and T is the symbol period. Generally, propagation in a radio channel is characterized by a very slowly varying time delay with respect to the channel gains affecting the signal amplitude and phase. Hence, it is reasonable to assume that the time delay remains constant and 978-1-5386-3531-5/17/$31.00 c 2017 IEEE the channel gain varies from one measurement to another. The received signal can thus be written as:

r(t) = α(t)s(t -τ ) + n(t), (2)
where τ is an unknown delay introduced by the channel, α(t) is a complex gaussian process representing the channel fading and n(t) is an additive white Gaussian noise (AWGN) of zero mean and variance σ 2 . The Nyquist pulse g(t) obtained from h(t), with g(t) = h(t) ⊗ h * (-t), satisfies the first Nyquist criterion (g(nT ) = 0 for any integer n = 0).

Let us consider the N -dimensional truncated vectors

r = [r 1 , • • • , r N ] T , s(τ ) = [s 1 (τ ), • • • , s N (τ )] T , α = [α 1 , • • • , α N ] T and n = [n 1 , • • • , n N ] T ,
representing the sampled version of r(t), s(t -τ ) and n(t), respectively taken at the rate F/T . α and τ are assumed to be unknown, respectively random and deterministic parameters. They can be gathered in the following unknown parameter vector:

µ = [α T , τ ] T . (3)
Let μ(r) be an unbiased estimator of the hybrid parameter µ = [α, τ ] T based on the observations r. The HCRB can be derived in an on-line context where the estimate of µ only depends on the current and the previous observations or in an off-line context when the whole observation block is received. This bound is defined as [START_REF] Trees | Detection, Estimation, and Modulation Theory: Radar-Sonar Signal Processing and Gaussian Signals in Noise[END_REF]:

E r,α|τ (μ(r) -µ) (μ(r) -µ) H ≥ HCRB(µ), (4)
where A ≥ B means that A -B is a positive semi definite matrix.

III. HCRB ANALYTICAL DERIVATION

From [START_REF] Rockah | Array shape calibration using sources in unknown locations-part i: Far-field sources[END_REF] and [START_REF] Hijazi | Bayesian Cramér-Rao Bound for ofdm rapidly time-varying channel complex gains estimation[END_REF], the HCRB can be decomposed into four block sub-matrices:

HCRB(µ) = HCRB 11 (α) hcrb 12 (α, τ ) hcrb 21 (α, τ ) hcrb 22 (τ ) . (5)
The HCRB is the inverse of the following hybrid information matrix (HIM) [START_REF] Rockah | Array shape calibration using sources in unknown locations-part i: Far-field sources[END_REF]:

H(µ) = E α [F (α, τ )] + E α|τ -△ µ µ log (P (α|τ )) , (6)
where P (α|τ ) is the a priori knowledge on α and F (α, τ ) is the Fisher information matrix (FIM) defined as:

F (α, τ ) = E r|α,τ -△ µ µ log (P (r|α, τ )) . (7)
P (r|α, τ ) is the conditional pdf of r given α and τ . Since α and τ are independent parameters, then P (α|τ ) = P (α) and thus:

E α|τ -△ µ µ log (P (α|τ )) = E α [-△ α α log (P (α))] 0 N,1 0 1,N 0 . (8)
The channel fading follows a Rayleigh distribution, therefore, α is a complex gaussian random variable of zero mean and covariance matrix R α . The a priori information on α is then given by:

P (α) = 1 |πR α | exp -α H R α -1 α . ( 9
)
By evaluating the second derivative of log (P (α)) with respect to α and making the expectation, we obtain:

E α [-△ α α log (P (α))] = R α -1 . (10)
The covariance matrix can be computed using the following expression [START_REF] Stüber | Principles of Mobile Communication[END_REF]:

[R α ] n,m = σ 2 α J 0 (2πf d T (n -m)) (11)
where σ α is the standard deviation of the channel fading, f d is the channel Doppler frequency and J 0 (.) is the Bessel function of the first kind.

The FIM can be decomposed into the following block matrices

F (α, τ ) = F 11 (α) f 12 (α, τ ) f 21 (α, τ ) f 22 (τ ) , (12)
where:

F 11 (α) = E r|α,τ [-△ α α log (P (r|α, τ ))] , (13) f 12 (α, τ ) = E r|α,τ [-△ τ α log (P (r|α, τ ))] , (14) f 21 (α, τ ) = E r|α,τ [-△ α τ log (P (r|α, τ ))] , (15) f 22 (τ ) = E r|α,τ [-△ τ τ log (P (r|α, τ ))] . (16
) In order to obtain the FIM, we need to first derive the likelihood of r expressed as, P (r|α, τ ).

A. Observation Likelihood Function

Based on the whiteness of noise samples, the likelihood function of r given α and τ is:

P (r|α, τ ) = N k=1 P (r k |α k , τ ), (17)
where:

P (r k |α k , τ ) = 1 2πσ 2 exp - |r k -α k s k (τ )| 2 2σ 2 . ( 18
)
Then ( 17) leads to:

P (r|α, τ ) = 1 2πσ 2 N exp - N k=1 |r k -α k s k (τ )| 2 2σ 2 .
(19) Making N tend to infinity, the observation likelihood function can be written in the continuous format:

P (r|α, τ )= 1 2πσ 2 exp -T0 |r(t) -α(t)s(t -τ )| 2 dt 2σ 2 = 1 2πσ 2 exp -T0 |r(t)| 2 + |α(t)s(t -τ )| 2 dt 2σ 2 + T0 ℜ {α(t) * s(t -τ ) * r(t)} dt σ 2 , ( 20
)
where T 0 is the observation period. We note that |r(t)| 2 does not depend on τ . According to [START_REF] Proakis | Digital Communications[END_REF], [START_REF] Mengali | Synchronization Techniques for Digital Receivers[END_REF], when T 0 is large enough, the impact of the time delay on the integral of the second squared term in [START_REF] Hijazi | Analytical Analysis of Bayesian Cramer-Rao Bound for Dynamical Rayleigh Channel Complex Gains Estimation in OFDM System[END_REF] can be neglected and thus:

T0 |α(t)s(t -τ )| 2 dt ≃ T0 |α(t)s(t)| 2 dt = N k=1 |α k | 2 |a k | 2 . ( 21
)
Considering:

x k (u) = y k (u) + v k (u), (22)
y k (u) = α k i aig ((k -i)T -(τ -u)) , (23)
v k (u) = T 0 h * (t -kT -u)n(t)dt, (24)
then after some mathematical derivations, the observation likelihood becomes:

P (r|α, τ ) = C exp - N k=1 |α k a k | 2 -2 N k=1 ℜ {α * k a * k x k (τ )} 2σ 2 , ( 25
)
where C is a constant term with respect to τ . x k (τ ) is seen as the matched filter output of the received signal taken at instant kT + τ (time delay perfectly recovered) , y k (τ ) is the useful part and v k (τ ) is a colored gaussian noise of zero mean and variance σ 2 . The symbols a k which appear in the updating term can be known if some pilot symbols are sent within the data frame. In this case, the TED is operating in a DA mode. In order to enhance the spectral efficiency, these pilot symbols can be limited to a short preamble or even omitted and a k are then estimated at the receiver. The TED is then operating in a Decision Directed (DD) mode. To do so, we usually implement a NDA mode in which a hard estimation of the symbols a k is used. Nevertheless, this technique performance degrades rapidly at low SNR, where modern systems are constrained to work.

Since we obtained the expression of the observation likelihood we can now compute the FIM following ( 13)- [START_REF] Yang | Bayesian and hybrid Cramer-Rao bounds for QAM dynamical phase estimation[END_REF].

B. Derivation of the Fisher information matrix in the DA mode

Based on ( 13)-( 16) and ( 25), we get:

F 11 (α) = diag(J 1 , • • • , J N ), (26)
where:

J k = |a k | 2 σ 2 , (27)
f 22 (τ ) = ρg(0), (28)
and

ρ = N k=1 |α k | 2 |a k | 2 σ 2
is the received SNR. We also obtain f 12 (α, τ ) = f 21 (α, τ ) T = 0N,1. This is due to the fact that α and τ are two independent parameters. In order to compute the analytical expression of the HIM given by ( 6), we further need to average the obtained FIM components over the vector α. Given that only f 22 (τ ) depends on α, then only the FIM related to the time delay τ has to be averaged over α. Averaging (28) over α is equivalent to averaging it over the SNR. We finally obtain:

Eα f 22 (τ ) = ρg(0), ( 29
)
where ρ is the SNR averaged over the fading channel.

C. Derivation of the Fisher information matrix in the NDA mode

In this paragraph, we provide the expression of the FIM components in the NDA mode for BPSK and QPSK modulated signals. In comparison with the DA mode, only the expression of f 22 (τ ) changes. For BPSK and QPSK signals, [START_REF] Mengali | Synchronization Techniques for Digital Receivers[END_REF]. The derivation of f 22 (τ ) for BPSK and QPSK modulations when the transmitted symbols are unknown at the receiver are given in [START_REF] Masmoudi | Closed-form expressions for the exact Cramer-Rao bounds of timing recovery estimators from BPSK, MSK and square-QAM transmissions[END_REF] for an AWGN channel by:

|a k | 2 = 1 in
f 22 (τ )=4ρ 1 - 1 √ 2π e -ρ β(ρ) ρ N m=1 N n=1 ġ((m -n)T ) - N 2 g(0) -ρ N m=1 N n=1 ġ((m -n)T ) , (30)
for a BPSK modulation where:

β(ρ) = +∞ -∞ e -x 2 2 cosh √ 2ρx dx, (31)
and:

f 22 (τ )= 2ρ 2 N m=1 N n=1 ġ2 ((m -n)T ) -2N ρg(0) × 1 √ 2π +∞ -∞ g 2 ρ (x) Gρ(x) e -x 2 2 dx - ρ π +∞ -∞ xgρ(x) 2 e -x 2 2 dx × N m=1 N n=1 ġ2 ((m -n)T ), (32)
for a QPSK modulation where:

g ρ (x) = 2e -2ρ sinh (2 √ ρx) , (33)
and

G ρ (x) = 2e -2ρ cosh (2 √ ρx) . (34)
Having a Rayleigh fading channel, then the SNR ρ is exponentially distributed and the pdf of ρ is given by:

P (ρ) = 1 ρ exp - ρ ρ , for ρ ≥ 0. ( 35
)
In order to obtain the final expression of the HIM, we need to average the analytical expressions of f 22 (τ ) for BPSK and QPSK modulations over the pdf of ρ for ρ ≥ 0. Since it is hard to obtain analytically a closed form expression of the integrals given by ( 30) and (32) and to make the averaging of the obtained result over the SNR distribution, a numerical integration can be used to obtain E α f 22 (τ ) . We note that the integrals in (30) and (32) decrease rapidly in terms of x. Thus, the integrand functions can be approximated by a finite Riemann integration over an interval [-D, +D] instead of ] -∞, +∞[ with an integration step δ. The same technique can be used to make the averaging over the SNR on the interval

[0, D] instead of [0, +∞[.

D. Inversion of the HIM

Based on (6) the HIM is given by:

H(µ) = Eα F 11 (α) + Rα -1 0 N,1 0 1,N Eα f 22 (τ ) (36)
We note that H(µ) is a block diagonal matrix thus using the block matrices inversion formula [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF], the HCRB is then given by:

HCRB(µ)= Eα F 11 (α) + Rα -1 -1 0 N,1 0 1,N Eα f 22 (τ ) -1 = HCRBα 0 N,1 0 1,N hcrbτ (37)
IV. SIMULATION RESULTS

In this section, we show simulation results of the HCRB for the channel coefficients and the time delay estimation. The diagonal elements of the inverse of HCRB(µ) give the HCRB with respect to α and τ . The results are provided for BPSK and QPSK modulated signals passed through a root raised cosine filter with a roll-off coefficient equal to 0.3. The up-sampling factor is equal to 8.

Figure 1 shows the HCRBs related to the channel and the time delay estimation versus the observation block length. HCRB α is given in an on-line context, where only the current and the previous samples are taken into account in the channel estimation problem. hcrb τ is computed for both DA and NDA modes and for BPSK and QPSK modulated signals. We have also fixed f d T = 0.001, which means that the channel varies slowly within the observation block. The online HCRB α is the same for the different modulations and the transmission modes since it only depends on the channel correlation matrix and the observation noise. We note that, like could be expected, the DA outperforms the NDA mode and that the BPSK modulation leads to lower time synchronization errors than the QPSK modulation. In Figure 2, we evaluate the on-line and the off-line HCRB α for various observation block lengths and for f d T = 0.001. Figure 3 displays, for a SNR= 5dB, the on-line and the off-line HCRB α for f d T = 0.1 (fast time varying channel) and f d T = 0.001 (almost invariant channel). It is shown that we can reach better estimation performance with an off-line technique, since we can use the whole observation block. The performances at the center of the symbols block are better than that of the borders for an off-line technique because the current, past and future observations are used at the same time for the channel gain estimation. However, at the borders we can only use the current and, according to the end or the beginning of the block, either the past or the future observation. It is also shown that, at the beginning, the on-line bound decreases, therefore the estimation performance is improved since the estimator uses more and more observations. After the reception of some signal samples, the bound converges to an asymptote and this is due to the observation noise. For instance, the asymptote is reached faster for f d T = 0.1 than for f d T = 0.001. In Figure 4, we display the evolution of the HCRB α with the SNR in the on-line and off-line scenarios and for different Doppler frequency values. It is shown that HCRB α decreases with the SNR and increases with fast varying channels. In Figures 5 and6, we evaluate, for respectively BPSK and QPSK modulations, the evolution of hcrb τ with the SNR for f d T = 0.1 and f d T = 0.001 and for various transmission modes (DA and NDA). The number of samples received within the signal block is equal to 40. As it could be expected, hcrb τ decreases as a function of SNR. Similar results to the previous cases are obtained by changing the SNR values.

V. CONCLUSION

In this contribution, we derive the analytical expression of the HCRB for time delay recovery under a random Rayleigh fading channel in the DA and the NDA modes. Simulation results have shown that we can reach better channel estimation using an off-line technique, where all the observation block is used. It is also shown that we can achieve a better time synchronization performance using the DA mode and with a slowly time varying channel.

Fig. 1 :

 1 Fig. 1: hcrb τ and HCRB α versus Observation block length (N ), f d T = 0.001, SNR=5dB.

Fig. 2 :Fig. 3 :

 23 Fig. 2: On-line and Off-line HCRB α versus Observation block length (N ), SNR=5dB.

Fig. 4 :

 4 Fig. 4: On-line and Off-line HCRB α versus SNR, N = 40.

Fig. 5 :

 5 Fig. 5: hcrb τ versus SNR, Block length= 40.

Fig. 6 :

 6 Fig. 6: hcrb τ versus SNR, Block length= 40.