G Farag

P Boivin

P Sagaut

Interaction of 2D spots with a heat releasing/absorbing shock wave: Linear Interaction Approximation results

The canonical interaction between a two-dimensional weak Gaussian disturbance (entropy spot, density spot, weak vortex) with an exothermic/endothermic planar shock wave is studied via the Linear Interaction Approximation. To this end, a unified framework based on an extended Kovasznay decomposition that simultaneously accounts for non-acoustic density disturbances along with a poloidal-toroidal splitting of the vorticity mode and for heat-release is proposed. An extended version of Chu's definition for the energy of disturbances in compressible flows encompassing multi-component mixtures of gases is also proposed. This new definition precludes spurious non-normal phenomena when computing the total energy of extended Kovasznay modes. Detailed results are provided for three cases, along with fully general expressions for mixed solutions that combine incoming vortical, entropy and density disturbances.

Introduction

The propagation of a hydrodynamic shock wave across an heterogeneous medium is a very important topic in many fields of application, e.g. aerospace engineering, nuclear engineering but also astrophysics. Such an interaction is known to emit a complex field, which is a mixture of acoustic, entropy and vortical waves according to Kovasznay's decomposition (see [START_REF] Kovasznay | Turbulence in supersonic flow[END_REF]; [START_REF] Chu | Non-linear interactions in a viscous heatconducting compressible gas[END_REF]; [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF]). In the limit of small disturbances, the emitted field can be accurately predicted considering a linearized theory, namely the Linear Interaction Approximation (LIA), see [START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF] for an exhaustive discussion. This approximation is relevant in the wrinkled shock regime, in which the shock front corrugation by upstream disturbances is small enough to leave its topology unchanged, so that it can be decomposed as a linear sum of sinusoidal contributions. Several semi-empirical criteria of validity of LIA have been proposed on the ground of Direct Numerical Simulation results. In the case of a turbulent upstream flow, [START_REF] Lee | Direct numerical simulation of isotropic turbulence interacting with a weak shock wave[END_REF] proposed

M 2 t < 0.1(M 2 1 -1) (1.1)
where M t and M 1 are the upstream turbulent and mean Mach numbers, respectively. This criterion was later refined using DNS with higher resolution by [START_REF] Ryu | Turbulence structure behind the shock in canonical shock-vortical turbulence interaction[END_REF], yielding

M t2 0.1M 2 (1.2)
with M t2 and M 2 the downstream (LIA-predicted) turbulent Mach number and the downstream mean flow based Mach number, respectively. In the laminar case of the interaction between an entropy spot and a normal shock, [START_REF] Fabre | Linear interaction of a cylindrical entropy spot with a shock[END_REF]) reported an excellent agreement within 1% error up to M 1 = 4 for disturbances with relative amplitude lesser or equal to 0.01. This theory was pioneered in the 1950s by Ribner (1954aRibner ( ,b, 1959)); [START_REF] Moore | Unsteady oblique interaction of a shock wave with a plane disturbance . Narita, Shinji 1973 The radiative energy loss from the shock front[END_REF] and is still under development. The most complete formulation of the normal mode analysis for canonical interaction was given by [START_REF] Fabre | Linear interaction of a cylindrical entropy spot with a shock[END_REF], which was further extended to the case of the non-reacting binary mixture of perfect gas [START_REF] Griffond | Linear interaction analysis applied to a mixture of two perfect gases[END_REF][START_REF] Griffond | A turbulent mixing reynolds stress model fitted to match linear interaction analysis predictions[END_REF] and to rarefaction waves (Griffond & Soulard 2012). Following this approach, wave vectors of emitted waves are obtained analytically thanks to the dispersion relation stemming from the linearized Euler equations, while wave amplitudes are solution of a linear system. A deeper physical insight is obtained by grouping upstream disturbances according to the Kovasznay normal-mode decomposition of small compressible fluctuations into acoustic, vorticity and entropy mode. This decomposition has been extended by splitting the vorticity mode as the sum of a poloidal and a toroidal components (Griffond & Soulard 2012), and also considering a binary mixture of perfect gas [START_REF] Griffond | Linear interaction analysis applied to a mixture of two perfect gases[END_REF][START_REF] Griffond | Linear interaction analysis for richtmyer-meshkov instability at low atwood numbers[END_REF]. Several cases have been succesfully investigated using LIA, among which the case of an upstream entropy spot [START_REF] Fabre | Linear interaction of a cylindrical entropy spot with a shock[END_REF], upstream vortical isotropic turbulent field [START_REF] Lee | Direct numerical simulation of isotropic turbulence interacting with a weak shock wave[END_REF][START_REF] Lee | Interaction of isotropic turbulence with shock waves: effect of shock strength[END_REF][START_REF] Quadros | Turbulent energy flux generated by shock/homogeneous-turbulence interaction[END_REF], upstream isotropic acoustic turbulent field [START_REF] Mahesh | The interaction of an isotropic field of acoustic waves with a shock wave[END_REF], upstream isotropic mixed vortical-entropy turbulent field [START_REF] Mahesh | The influence of entropy fluctuations on the interaction of turbulence with a shock wave[END_REF]).

An alternative complete analytical treatment of the linearized problem based on the Laplace transform has been developed by Wouchuk, Huete and coworkers in a series of papers (e.g. [START_REF] Wouchuk | Analytical linear theory for the interaction of a planar shock wave with an isotropic turbulent vorticity field[END_REF][START_REF] De Lira | Turbulence generation by a shock wave interacting with a random density inhomogeneity field[END_REF]Huete et al. 2012aHuete et al. ,b, 2013)). Here, a telegraphist equation is obtained for each type of incident wave whose analytical solution gives the amplitude of emitted disturbances. This approach has not been explicitly recast into Kovasznay framework up to now, but acoustic and vortical upstream fluctuations have been considered in a series of papers, along with density fluctuations. The analysis has been recently extended to the case of thin detonation waves [START_REF] Huete | Theory of interactions of thin strong detonations with turbulent gases[END_REF][START_REF] Huete | Linear theory for the interaction of small-scale turbulence with overdriven detonations[END_REF], which are described as shock wave associated to a heat release phenomenon. That approach has also been applied to many cases, e.g. incident isotropic adiabatic turbulence [START_REF] Wouchuk | Analytical linear theory for the interaction of a planar shock wave with an isotropic turbulent vorticity field[END_REF], pure incident acoustic turbulence (Huete et al. 2012b), pure incident isotropic density fluctuations including the re-shock problem (Huete et al. 2012a).

Selected studies carried out within these two general frameworks are listed in Table 1 in an attempted summary, sorting the studies referred to in the two previous paragraphs according to the perturbation modes considered, the possibility to account for heat releasing/absorbing shock, as well as the upstream perturbations and the approach followed. It is worth noting that in the case of an upstream turbulent field, LIA can be rewritten in terms of turbulent fluxes, leading to a linear problem for the jump of these quantities across the shock. These relations can be used to derive RANS models well suited for the simulation of the shock-turbulence interaction [START_REF] Sinha | Modeling shock unsteadiness in shock/turbulence interaction[END_REF][START_REF] Griffond | A turbulent mixing reynolds stress model fitted to match linear interaction analysis predictions[END_REF][START_REF] Soulard | Pseudocompressible approximation and statistical turbulence modeling: Application to shock tube flows[END_REF][START_REF] Sinha | Evolution of enstrophy in shock/homogeneous turbulence interaction[END_REF][START_REF] Quadros | Turbulent energy flux generated by shock/homogeneous-turbulence interaction[END_REF].

The goal of the present paper is three-fold. First, it aims at providing a complete, unified formulation of the normal-mode-based LIA approach that encompasses all previous developments, namely binary mixture of perfect gas interacting with a nonadiabatic shock wave considering the poloidal/toroidal splitting of vorticity. The various extensions mentioned above have not been gathered into a single unified framework up to now. In particular, accounting for the non-adiabatic character of a shock wave simultaneously with these extensions has not been done up to now, although it was (turbulent) and Spot refer to the nature of the upstream field. The approach followed is also indicated as L/O, referring respectively studies articles based/not based on the Laplace transform. carried out in the case of density fluctuations through detonations [START_REF] Huete | Theory of interactions of thin strong detonations with turbulent gases[END_REF][START_REF] Huete | Linear theory for the interaction of small-scale turbulence with overdriven detonations[END_REF]. Heat-release/absorption will be described as a punctual source/sink at the shock, to encompass thin reactive shock waves, shock-induced condensation or radiative loss (see e.g. Zel 'Dovich & Raizer 2012). In this general formulation, all types of upstream disturbances will be considered within an extended Kovasznay decomposition framework.

The second goal of the paper is to extend Chu's definition for disturbance energy [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (part i)[END_REF]) to a multi-component fluid: a physically relevant and mathematically consistent definition well-suited for small perturbations definition of the disturbance energy is of primary importance to analyze the effect of the interaction with the shock wave, and is therefore a prerequisite to the next paper's goal.

The last goal of the present paper is to analyze the interaction of a Gaussian perturbation spot with a shock wave in the presence of phenomena mentioned above. Three different cases are investigated: a density spot, an entropy spot and a vorticity spot (i.e. a weak vortex). It is worth noting that the case of upstream density heterogeneities has been considered in the case of non-reactive shock waves and thin strong detonations by [START_REF] Huete | Theory of interactions of thin strong detonations with turbulent gases[END_REF]. Such a simple configuration can be considered as an idealized model of the interaction of a shock wave with a two-phase heterogeneity (bubble, droplet) with small density ratio. To the knowledge of the authors, such general cases have never been considered in the open literature up to now. Using the three elementary cases considered in the present papers, an infinite number of cases can be derived by linear combination of the LIA results. As an example, the interaction between a shock wave and a cold weak vortex is obtained in a straightforward way by linear combination of the solutions related to the isentropic vortex case and a cold entropy spot. Multiple spot solutions can also be found in the same way, introducing a space-time shift in the solution associated to each spot. The optimal combination of these elementary spots to minimize the radiated noise is investigated in the present paper.

The paper is organized as follows. The basic physical model and associated governing equations are displayed in Section 2. The decomposition of both upstream and downstream fields according to the present extended Kovasznay modal decomposition is then presented in Section 3. The extended definition of disturbance energy and its relation to the energy of Kovasznay modes are discussed in Section 4. Then the proposed general formulation of the normal-mode-based LIA approach is discussed in Section 5. The application to the interaction of a heat releasing/absorbing shock wave with a variety of Gaussian spots (for density, entropy and vorticity fluctuations) is then addressed in Section 6, with most of the technical details regarding the treatment of 2D Gaussian spots given in Appendix. Conclusions are drawn in Section 7.

Physical model

The physical model addressed in the present paper is related to the case of 2D canonical shock/disturbance interaction in a binary mixture of perfect gas, in the presence of heat release/absorption on the shock wave. Viscous effects are neglected. Upstream and downstream of the normal shock, the flow is governed by the Euler equations:

                     ∂ρ ∂t + div(ρu) = 0, ∂ρu ∂t + div ρu ⊗ u + pI = 0, ∂ρE ∂t + div ((ρE + p)u) = 0, ∂ρY ∂t + div(ρY u) = 0, (2.1)
where p, ρ, u and E denote the mixture pressure, density, velocity, and total energy; and Y is the mass fraction of the first component in the binary mixture (see e.g. [START_REF] Williams | Combustion Theory[END_REF].

The mixture equation of state for the binary mixture reads

p = ρ R W T , (2.2)
where R and W denote the perfect gas constant and the molar weight of the mixture, respectively. The classical relations for ideal gas mixtures yields the following relations between the component properties and the mixture properties:

1/W = Y /W a + (1 -Y )/W b , c v = Y c v,a + (1 -Y )c v,b , and c p = Y c p,a + (1 -Y )c p,b , γ a = c p,a c v,a , γ b = c p,b c v,b , γ = c p c v , A r t = r a -r b r , and A cv t = c va -c vb cv . (2.3)
W , c v , c p and γ denote respectively the mixture molecular weight, mass heat capacity at constant volume and constant pressure, the heat capacity ratio, as well as two Atwood numbers, to be used hereafter. Subscripts a and b denote the corresponding component thermodynamic properties in the binary mixture, one being inert and one possibly reactive. Note however that they do not intervene in the following, where indices exclusively serve as to identify the upstream and downstream states.

Considering the case of an 1D flow along the x axis and a normal shock wave and denoting (u x , u r , u φ ) the component of velocity in cylindrical coordinates (in the discontinuity reference frame, the x axis being taken normal to the planar shock wave), the upstream and downstream mean quantities (resp. subscripts 1,2) relate through the Hugoniot jump conditions for mass, momentum and energy:

ρ 1 u x1 = ρ 2 u x2 , p 1 + ρ 1 u 2 x1 = p 2 + ρ 2 u 2 x2 , h 1 + u 2 x1 2 = h 2 + u 2 x2 2 , (2.4)
with u r , u φ and Y a being conserved through the shock:

u r1 = u r1 , u φ1 = u φ2 , Y a,1 = Y a,2 .
The enthalpy h jump condition may be reformulated as

c p T 1 + u 2 x1 2 = c p T 2 + u 2 x2 2 -∆Q, (2.5)
where ∆Q accounts for heat release/heat absorption at the shock wave. ∆Q > 0 was considered by [START_REF] Huete | Theory of interactions of thin strong detonations with turbulent gases[END_REF] to model thin detonations, while ∆Q < 0 should be used to account for physical mechanisms restricted to a thin region downstream the shock front that act as an energy sink, e.g. radiative losses or condensation [START_REF] Zel'dovich | Physics of shock waves and high-temperature hydrodynamic phenomena[END_REF].

Note that, while ∆Q is here formulated as an independent parameter, a classical assumption for strong detonations (see, e.g. [START_REF] Williams | Combustion Theory[END_REF] , the heat absorption typically depends on the shock strength for endothermic processes (which typically ends when saturation is reached), as is the case in ionizing, nuclear dissociating shocks as those occurring in core collapsing supernovae [START_REF] Huete | The impact of vorticity waves on the shock dynamics in core-collapse supernovae[END_REF][START_REF] Abdikamalov | Turbulence generation by shock-acoustic-wave interaction in core-collapse supernovae[END_REF][START_REF] Huete | Response of nuclear-dissociating shocks to vorticity perturbations[END_REF], shock-induced condensation in vapor-liquid two-phase flow [START_REF] Zhao | Shock wave of vapor-liquid two-phase flow[END_REF] or cooling induced by radiative loss (Narita 1973).

Introducing the sound speeds on either side of the shock c 1 , c 2 in the jump conditions lead to the following relation between the upstream and downstream Mach Numbers, respectively M 1 and M 2 :

M 2 2 = 1 + γM 2 1 -(M 2 1 -1) √ 1 -β 1 + γM 2 1 + γ(M 2 1 -1) √ 1 -β , β = 2(γ 2 -1)M 2 1 q (1 -M 2 1 ) 2 , (2.6)
where the normalized heat coefficient has been introduced

q = ∆Q c 2 1 .
(2.7)

The compression factor m = ρ 2 /ρ 1 = u 1 /u 2 is obtained through

1 m = 1 + γM 2 1 (γ + 1)M 2 1 + 1 -M 2 1 (γ + 1)M 2 1 1 -β.
(2.8)

Note that c p and γ, appearing in the above relations are identical on both side of the shock thanks to the continuity of mass fraction Y , thereby considerably reducing the equations. A detailed account on the validity of this assumption has been provided by [START_REF] Griffond | Linear interaction analysis applied to a mixture of two perfect gases[END_REF]: the analysis is valid for small concentration fluctuations within binary mixtures with very different thermodynamic properties, or large concentration fluctuations within gases of similar thermodynamic properties. This translates, in practice, to the assumption holding when the reactive component mixture is sufficiently dilute in the inert one, as is often the case in air. When the assumption does not hold, the present study still present valuable benchmarks for numerical codes, in which thermodynamic properties may be artificially set to constants. All other classical relations for T 2 /T 1 , p 2 /p 1 , ... are formally identical to those of the classical normal shock case, M 2 and m being now given by the above formula.

The consistency constraint which ensures that both m and M 2 remain positive is q min < q < q max , (2.9)

where

q min = 1 1 -γ - M 2 1 2 , q max = (1 -M 2 1 ) 2 2(γ 2 -1)M 2 1 .
(2.10)

The consistent domain for heat-source/sink as a function of the upstream Mach number M 1 is illustrated in Fig. 1. Superimposed are contours of the downstream Mach number M 2 , as provided by (2.6). One recovers the physical behavior that the downstream flow is accelerated in the case q > 0 compared to the neutral shock case q = 0, while it is decelerated in the opposite case q < 0 , due to the balance between kinetic energy and internal energy. In the asymptotic limit q = q max , the system satisfies the so-called Chapman-Jouguet condition M 2 = 1 (see, e.g. [START_REF] Zeldovich | On the theory of the propagation of detonation in gaseous systems[END_REF]). The other limit, q = q min corresponds to an infinite mass compression ratio, impossible to sustain in practice. For this reason, the endothermic cases presented in Section 6 are presented for q = q min /2, translating to at most half the upstream kinetic energy being absorbed, leading to reasonable compression ratio and downstream Mach numbers (resp. m = 6.5 and M 2 = 0.33).

The Kovasznay modal decomposition for disturbances in a binary mixture of ideal gas

The Linear Interaction Approximation relies on a small disturbance hypothesis and the use of linearized equations to described fluctuation propagation on either side of the shock.

For each quantity (e.g. u), let us identify the fluctuation part (u ′ ) and the mean (ū) as

u = ū + u ′ , p = p + p ′ , . . .
and assume the fluctuation part is small (u ′ /ū ≪ 1), a classical assumption provided:

• Linearization of Y , for which Ȳ = 0 is acceptable, is valid [START_REF] Griffond | Linear interaction analysis applied to a mixture of two perfect gases[END_REF]. This is in practice related to the continuity of c p and γ discussed after Eq. (2.8).

• Similarly, the linearization for the normal shock velocity is questionable in the limit ū → 0, attainable when q → q min . To avoid this, the present study should not be carried out for M 2 < 0.25, or, alternatively, q < q min /2.

In the reference frame tied to the planar shock front the 2D perturbation field then satisfies

                         ∂ρ ′ ∂t + ū ∂ρ ′ ∂x + ρ ∂u ′ j ∂x j = 0, ∂u ′ i ∂t + ū ∂u ′ i ∂x + 1 ρ ∂p ′ ∂x i = 0, ∂Y ′ ∂t + ū ∂Y ′ ∂x = 0, ∂p ′ ∂t + ū ∂p ′ ∂x + γ p ∂u ′ j ∂x j = 0, (3.1)
which can be recast as a system of evolution equations for Kovasznay's physical modes:

                               ∂s ′ ∂t + ū ∂s ′ ∂x = 0, ∂ω ′ ∂t + ū ∂ω ′ ∂x = 0, ∂ω ′ ⊥ ∂t + ū ∂ω ′ ⊥ ∂x = 0, ∂ ∂t + ū ∂ ∂x 2 p ′ = c 2 ∇ 2 p ′ , ∂Y ′ ∂t + ū ∂Y ′ ∂x = 0, (3.2)
where ω ′ = ∇×u ′ denotes the fluctuating vorticity, and

ω ′ ⊥ = (ω ′ •n)n and ω ′ = ω ′ -ω ′ ⊥
are the shock-normal and the shock-parallel components of vorticity, respectively, with n the unit normal vector of the planar shock wave. The shock-normal and the shocktangential components correspond to the toroidal and poloidal components of the velocity field in the reference frame tied to the planar shock front, respectively. One recognizes the entropy mode, the toroidal and poloidal vorticity modes, the fast and slow acoustic modes and the concentration mode. It is worth noting that Kovasznay's modes correspond to the eigenmodes of the linearized propagation operator, which are orthonormal according to the inner product associated to Chu's definition of compressible disturbance energy.

k a k s k x y α a α s α x s (y, t)
Figure 2: Sketch of the configuration. The corrugated shock mean front position is at x = 0. The incident perturbation has wave vector k, at angle α with respect to the shock normal. Emitted waves may be acoustic waves, with wave vector k a , or non-acoustic ones, with wave vector k s .

Let us now introduce propagating plane wave disturbances of the general form

φ ′ = A i (k) exp [i(k • x -Ωt)].
(3.3)

Here, A i (k) denotes the amplitude of upstream Kovasznay mode of type i, with i = s, a, Y, v, t for entropy, acoustic, concentration and poloidal/toroidal vorticity mode, respectively. k is the perturbation wave vector, associated with pulsation Ω = ū1 k cos α, where α is the angle of the incident perturbation with respect to the shock, as illustrated in Fig. 2. The upstream fluctuating field can then be decomposed as follows

            τ ′ 1 /τ 1 u ′ x1 /ū 1 u ′ r1 /ū 1 u ′ φ1 /ū 1 p ′ 1 /γ p1 Y ′ 1 T ′ 1 / T1 s ′ 1 /C p1             = A i (k)e i(k•x-Ωt)             δ is -δ ia + δ iY A r t δ iv sin α + δ ia cos α M1 -δ iv cos α + δ ia sin α M1 δ it δ ia δ iY δ is + (γ -1)δ ia δ is             , (3.4)
where δ ij is the Kronecker symbol, and τ = 1/ρ is the specific volume. Now introducing the transfer function Z ij between upstream Kovasznay mode of type i and downstream Kovasznay mode of type j, the emitted fluctuating field downstream the shock is given by:

            τ ′ 2 /τ 2 u ′ x2 /ū 2 u ′ r2 /ū 2 u ′ φ2 /ū 2 p ′ 2 /γ p2 Y ′ 2 T ′ 2 / T2 s ′ 2 /C p2             = A i (k)e -ka ηx e i(ka•x-Ωt)             -Z ia Z ia (cos α a + iη)/(M 2 ζ) Z ia sin α a /(M 2 ζ) 0 Z ia 0 (γ -1)Z ia 0             + A i (k)e i(ks•x-Ωt)             Z is + Z iY A r t Z iv sin α s -Z iv cos α s Z it 0 Z iY Z is Z is             , (3.5)
where

ζ = 1 -η 2 + 2iη cos α a . (3.6)
Acoustic and non-acoustic emitted fluctuations are separated into two contributions in Eq. (3.5), as they correspond to different wave vectors, resp. k a (possibly associated to attenuation η) and k s . These wave vectors are explicited hereafter. The transfer function Z ij coefficients are explicitly given in Section 5.

Emitted acoustic and non-acoustic wave vectors

Evaluation of the wave vectors k a , k s and the associated angles α a , α s and attenuation η is classical (see [START_REF] Fabre | Linear interaction of a cylindrical entropy spot with a shock[END_REF][START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF], but is nonetheless recalled for the sake of completness.

The effect of α being different whether the incident perturbation is acoustic (i = a) or non-acoustic (i = a), it is convenient to introduce the modified incident angle β as

β = δ ia α ′ + (1 -δ ia )α, (3.7)
where α ′ is defined as

cot α ′ = cot α + 1 M 1 sin α . (3.8)
Wave vectors and angles are then related through the relation:

k a,s k = sin β sin α a,s , (3.9)
valid for both acoustic and non-acoustic emitted waves.

The emitted non-acoustic wave vector angle simply reads

cot α s = m cot β, (3.10)
where m is the compression factor (2.8).

Obtaining the emitted acoustic wave vector k a and associated attenuation η is not as straightforward. If the incident perturbation is non-acoustic (i = a), a singularity appears for α = β = α c , where

cot α c = 1 -M 2 2 mM 2 , (3.11)
for which the emitted acoustic wave vector corresponds to the critical emission angle cos α c a = -M 2 .

(3.12)

If α < α c , the emitted wave vector angle reads

cot α a cot α c a = cot β cot α c - 1 M 2 cot β cot α c 2 -1 and η = 0, (3.13) else if α > α c : cot α a cot α c a = cot β cot α c and η = | cot α c a sin α a | M 2 - 1 M 2 1 - cot β cot α c 2 , (3.14)
In the particular case where the incident perturbation is acoustic (i = a), β = α ′ , and two critical values for the incident angle α are found:

α ± c + 1 M 1 sin α ± c = ∓ 1 -M 2 2 mM 2 , (3.15)
corresponding to fast and slow propagation regimes, separated by incident angle α M such as

cos α M = - 1 M 1 . (3.16)
For acoustic incident perturbations (β = α ′ ) Eqs. (3.13) and (3.14) remain valid, now defining four regimes: a propagating and non-propagating regime for each of the fast and slow modes.

The global procedure for the determination of emitted wave vectors as well as associated attenuation is summarized in Table 2 andthe

             ∀α, αs from Eq. (3.10) with β = α ′ 0 < α < α - c , (αa, η) from (3.13) with β = α ′ α - c < α < αM , (αa, η) from (3.14) with β = α ′ αM < α < α + c , (αa, η) from (3.14) with β = α ′ α + c < α < π, (αa, η) from (3.13) with β = α ′
Table 2: Computation of the emitted acoustic and non acoustic wave vectors through the corresponding angles α a and α s , for non-acoustic and acoustic incident perturbation. Also included is the determination of attenuation η for emitted acoustic waves. angle is illustrated in Fig. 4.

Extension of Chu's definition for disturbance energy to multicomponent gas

An important issue is the derivation of a physically relevant and mathematically consistent definition of the energy of the disturbances in compressible flows. Chu's definition [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (part i)[END_REF] for the disturbance energy around a base flow has the advantage to define an inner product, with respect to which the linearized Euler equations about a uniform base flow are self-adjoint, and Kovasznay modes correspond to orthogonal eigenmodes of the linearized operator. The orthogonality of eigenmodes prevents spurious non-normality-induced phenomenon in the computation of the energy of the fluctuating field [START_REF] George | On chu's disturbance energy[END_REF][START_REF] Sagaut | Homogeneous turbulence dynamics[END_REF] As a matter of fact, the use of a non-normal basis may lead to unphysical growth of the energy of the system because of the contributions of non-zero cross-products of basis vectors. Therefore, one can split the total energy as the sum E tot = i E i , with i = v, a, s for the vorticity mode, the acoustic mode and the entropy mode, respectively.

Since the present work deals with multi-component gas, the original Chu's definition is extended in the present section. A first step consists of finding an expression of the linearized Euler equations that will lead to orthogonal eigenvectors. This is the case when the matrix associated to the linearized problem is symmetric. To this end, an adequate choice of physical unknowns must be done. Noticing that the set (ρ, u, v, T, Y ) leads to a non-symmetric matrix and non-orthogonal eigenvectors, we choose here to write the linearized problem using (ρ a , ρ b , u, v, T ) :

                                   ∂ρ ′ a ∂t + ū ∂ρ ′ a ∂x + ρa ∂u ′ j ∂x j = 0, ∂ρ ′ b ∂t + ū ∂ρ ′ b ∂x + ρb ∂u ′ j ∂x j = 0, ∂u ′ i ∂t + ū ∂u ′ i ∂x + 1 ρ ∂p ′ ∂x i = 0, ∂T ′ ∂t + ū ∂T ′ ∂x + p ρC v ∂u ′ j ∂x j = 0, p ′ p - T ′ T -[1 + A r t (1 -Y 0 )] ρ ′ a ρ 0 -(1 -A r t Y 0 ) ρ ′ b ρ 0 = 0, (4.1)
where the last line is related to the linearized equation of state, with

ρ a = ρY ρ b = ρ(1 -Y ), (4.2) ρ ′ a ρ 0 = ρ ′ ρ 0 Y 0 + Y ′ , (4.3) ρ ′ b ρ 0 = ρ ′ ρ 0 (1 -Y 0 ) -Y ′ , (4.4)
Now introducing the vector of normalized variables X = ( ρa , ρb , ũ, ṽ, T ) T where

ũ = u ′ c 0 ; ṽ = v ′ c 0 ; T = T ′ T 0 γ(γ -1) ; ρa = ρ ′ a ρ 0 γY 0 1 + A r t (1 -Y 0 ) ; ρb = ρ ′ b ρ 0 γ(1 -Y 0 ) 1 -A r t Y 0 (4.5)
and considering propagating plane wave disturbances, the linearized problem (4.1) can be rewritten in the following compact form

dX dt = M X, (4.6)
where the linearized operator matrix is given by

M =          -ik x u 0 0 -ikxc0 √ γ K 1 - ikyc0 √ γ K 1 0 0 -ik x u 0 -ikxc0 √ γ K 2 - ikyc0 √ γ K 2 0 -ikxc0 √ γ K 1 -ikxc0 √ γ K 2 -ik x u 0 0 -ikxc0 √ γ-1 √ γ - ikyc0 √ γ K 1 - ikyc0 √ γ K 2 0 -ik x u 0 - ikyc0 √ γ-1 √ γ 0 0 -ikxc0 √ γ-1 √ γ - ikyc0 √ γ-1 √ γ -ik x u 0          (4.7)
where the two positive parameters K 1 and K 2 are defined as

K 1 = Y 0 [1 + A r t (1 -Y 0 )], K 2 = (1 -Y 0 )(1 -A r t Y 0 ), (4.8)
The five eigenvalues are

-ik x u 0 , -ik x u 0 , -ik x u 0 , -i (k x u 0 ∓ kc 0 ) , (4.9)
which correspond to the normalized propagation speeds of (from the left to the right) the entropy mode, the vorticity mode, the concentration mode and the fast and slow acoustic modes. The associated set of orthogonal eigenvectors is

X s =         γ-1 γ+K 2 1 -1 0 0 0 -K1 √ γ+K 2 1 -1         , X v =       0 0 ky k -kx k 0       , X a ± =         K1 √ 2γ K2 √ 2γ ±kx k √ 2 ±ky k √ 2 γ-1 2γ         , (4.10) X Y =          -K1K2 √ (γ+K 2 1 -1)γ γ+K 2 1 -1 γ 0 0 -K 2 γ-1 (γ+K 2 1 -1)γ          . (4.11)
All possible solutions of the linearized problem can be expressed as a linear combination of the eigenvectors:

X(t) = i=s,v,a ± ,Y C i (t)X i .
Therefore a local definition of the total energy E(t) of the disturbance is given by the square of L 2 norm of X(t). Thanks to the orthogonality property, one has

X(t) 2 = X(t) • X(t) = i=s,v,a ± ,Y C 2 i (t) X i 2
, which appears as the sum of the energy of each mode. The associated energy in a volume V is obtained in a straightforward way as:

E tot (t) = γp 0 2 V K 2 1 ρ ′ 2 a γρ 2 0 Y 2 0 + K 2 2 ρ ′ 2 b γρ 2 0 (1 -Y 0 ) 2 + u ′ i u ′ i c 2 0 + T ′ 2 γ(γ -1)T 2 0 dV , (4.12)
which can be rewritten as a function of u ′ i , p ′ , s ′ and Y ′ as follows:

E tot (t) = γp 0 2 V M 2 0 u ′ i u 0 2 + p ′ γp 0 2 + 1 γ K 2 1 Y 2 0 + K 2 2 (1 -Y 0 ) 2 + (A r t ) 2 γ -1 (Y ′ ) 2 (4.13) + 1 γ -1 s ′ C p 2 dV.
The original formula given by Chu for single-species fluids is recovered taking Y 0 = 1 (which leads to K 1 = 1, K 2 = 0) along with Y ′ = 0.

A general formulation of the normal-mode-based LIA

The shock jump relations for a normal planar shock wave with possible heat release/absorption and change in specific heats across the shock read

ρ1 (u ′ x1 - ∂x s ∂t ) + ū1 ρ ′ 1 = ρ2 (u ′ x2 - ∂x s ∂t ) + ū2 ρ ′ 2 ,
(5.1)

p ′ 1 + ρ ′ 1 ū2 1 + 2 ρ1 ū1 u ′ x1 = p ′ 2 + ρ ′ 2 ū2 2 + 2 ρ2 ū2 u ′ x2 , h ′ 1 + ū1 (u ′ x1 - ∂x s ∂t ) = h ′ 2 + ū2 (u ′ x2 - ∂x s ∂t ), ū1 ∂x s ∂y + u ′ r1 = ū2 ∂x s ∂y + u ′ r2 , u ′ φ1 = u ′ φ2 , Y ′ 1 = Y ′ 2 .
As in Eq. (3.4), all prime quantities (e.g. p ′ 1 ) correspond to the fluctuations around the average base flow (e.g. p1 ), and

x s = x s (y, t) = A x e i(k sin βy-Ωt ),
(5.2) denotes the shock displacement with respect to its equilibrium position, as depicted in Fig. 2. A x is the perturbation amplitude.

The jump relations (5.1) can be normalized as

- τ ′ 2 τ2 + u ′ x2 ū2 -i cos β(1 -m) A x = - τ ′ 1 τ1 + u ′ x1 ū1 , - τ ′ 2 τ2 + 2 u ′ x2 ū2 + 1 M 2 2 p ′ 2 γ p2 = m - τ ′ 1 τ1 + 2 u ′ x1 ū1 + 1 M 2 1 p ′ 1 γ p1 , (5.3) u ′ r2 ū2 + i sin β(1 -m)A x = m u ′ r1 ū1 , u ′ x2 ū2 + 1 M 2 2 + 1 (γ -1) M 2 2 p ′ 2 γ p2 + 1 (γ -1) M 2 2 τ ′ 2 τ2 + i cos β(1 -m) m A x = m 2 u ′ x1 ū1 + 1 M 2 1 + 1 (γ -1) M 2 1 p ′ 1 γ p1 + 1 (γ -1) M 2 1 τ ′ 1 τ1
,

+ 1 γ(γ -1) m 2 M 2 1 - 1 M 2 2 A Cv t -A r t , u ′ φ2 ū2 = m u ′ φ1 ū1 , Y ′ 1 = Y ′ 2 .
From the normalized shock relations, the transfer functions introduced in (3.5) can be expressed through the linear system (5.4) where the transfer function vector Z i contains the intensity of each emitted Kovasznay mode for a given incident mode i = Y, t, v, s, a, x

M Z i = B i ,
Z i = (Z iY , Z it , Z iv , Z is , Z ia , Z ix ) T .
(5.5)

The matrix M reads

M =           1 0 0 0 0 0 0 1 0 0 0 0 0 0 sin α s -1 1 + cos αa+iη M2ζ i(m -1) cos β 0 0 2 sin α s -1 M 2 2 +1 M 2 2 + 2 cos αa+iη M2ζ 0 0 0 -cos α s 0 sin αa M2ζ i(1 -m) sin β 0 0 sin α s 1 (γ-1)M 2 2 1 M 2 2 + cos αa+iη M2ζ im (1 -m) cos β           , (5.6)
and the right-hand term, dependent on the incident wave's nature:

B s =         0 0 -1 -m 0 m 2 (γ-1)M 2 1         , B v =         0 0 sin α 2m sin α -m cos α m 2 sin α         , B a =           0 0 1 + cos α M1 m M 2 1 +1 M 2 1 + 2 cos α M1 m sin α M1 m 2 1 M 2 1 + cos α M1           , B t =         0 m 0 0 0 0         , B Y =          1 0 0 (1 -m)A r t 0 m 2 M 2 1 -1 M 2 2 1 γ(γ-1) A Cv t + 1 γ A r t         
.

(5.7) From the above system, the transfer function vector can be deduced as

(Z iY , Z it , Z iv , Z is , Z ia , Z ix ) T = M -1 B i , (5.8)
where the inverse matrix is a block diagonal matrix of the same form as M . It can then be inferred that the toroidal mode is fully decoupled from the others

Z tj = B t , Z it = 0 for i = t.
(5.9)

A similar behavior is obtained for the concentration mode Y , when (A r t , A Cv t ) = (0, 0) and B Y comprises of a single non-zero component. For arbitrary values of (A r t , A Cv t ), however,

Z Y j = B Y , Z iY = 0 for i = Y,
(5.10) so that an upstream mass concentration perturbation can produce a combination of various modes downstream of the shock. Downstream, however, a mass fraction perturbation can only arise from an upstream mass fraction perturbation. These comments allow to consider a reduced number of Z ij terms in the following Figures.

The transfer functions obtained for acoustic, poloidal and entropy incident perturbations are plotted in Fig. 5 as functions of the incident angle α. The associated emitted wave vectors are found in Fig. 4 (α a for Z ai and α s for Z vi and Z si ).

Incident mass fraction perturbations can vary in nature depending on the value of Atwood's numbers (A r t , A Cv t ) defined earlier (2.3). The associated transfer function Z Y i are therefore provided separately, in Fig. 6, with associated emitted wave angle α s in Fig. 4. Note that B Y is linear in A r t and A Cv t , so that providing solutions Z Y i for the two base vectors (A r t , A Cv t ) = (0, 1) and (A r t , A Cv t ) = (1, 0) suffice to describe the transfer function for any (A r t , A Cv t ).

Interaction with Gaussian spots

This section is dedicated to the interaction between 2D Gaussian spots advected at the uniform speed U 1 in the shock-normal direction and a planar shock wave. The Gaussian spots are introduced as perturbations of the form

G ′ = ǫ e -r 2 , (6.1)
where r is the radial coordinate relative to the centre of the spot, and the Gaussian perturbation G ′ is successively set as three elemental perturbations

                 G ′ = τ ′ A r t τ
for the density spot,

G ′ = s ′ c p
for the entropy spot,

G ′ = ω ′ U
for the vorticity spot.

(6.2)

For each perturbation, the emitted flow will systematically be studied through comparisons of acoustic, entropy and vorticity fields. Note that, owing to the linear character of this study, it is straight-forward to combine these three elemental Gaussian perturbations into more complex ones, and obtain the emitted flow-field. This is illustrated in Fig. 7, which displays the vorticity field emitted from the combination of the three elemental spots presented hereafter.

In the following ǫ, appears as a mere scaling and is therefore set to 1. Typical results are shown for A r t = 2, A Cv t = 1, M 1 = 2 and γ = 1.4. To illustrate the effect of the heat-release, results are plotted for adiabatic (q = 0), endothermic (q = -2.25) and exothermic (q = 0.59). The numeric values for endothermic and exothermic shocks were chosen to be q min /2 and q max /2 at M 1 = 2.

Gaussian density spot

Let us now consider a density spot, e.g . G ′ = τ ′ A r t τ in (6.2), which can be considered as an idealized model for shock/dropplet interaction.

The choice of a positive ǫ in (6.2) corresponds to the definition of a heavy perturbation with respect to the upstream fluid, which can be interpreted as an ideal model for a droplet of heavy fluid. A negative value would correspond a pocket of light fluid. It is x x q > 0 q = 0 worth noting that pure density heterogeneities without acoustic perturbation, i.e. pure ρwaves, are obtained considering concentration fluctuations. The solution is then computed analytically thanks to the formulas given in the Appendix. The emitted fields of normalized entropy s ′ Cp and vorticity aω ′ ū are displayed in the first 4 plots of Figure . 8. Since the emitted patterns are advected at the constant speed U 2 , they are plotted in the reference frame associated to the perturbation centre, in which they are frozen thanks to the fact that diffusive effects are not taken into account in the present inviscid model. The presented patterns are related to the far field solution, i.e. intermediary solutions that are found at times at which the incoming fluctuation spot has not totally crossed the shock are not presented for the sake of brevity (but can be computed).

-7 -3 -3 2 2 2 10 0 2 4 -2 -1 0 1 2 -1 -0 .6 -0 .2 -0 .2 0. 2 0.2 0.2 0 .2 0. 6 0.6 2 4 -4 -2 0 -2 -1 0 1 2 x x q < 0 aω ′ ū -0 .3 -0 .1 -0 .1 0 .0 2 5 0.0 25 0.025 0 .0 2 5 0.2 0.2 1 1 0 2 4 -2 -1 0 1 2 -1 -0 .6 -0 .2 -0 .2 0. 2 0.2 0.2 0 .2 0. 6 0.6 2 4 -4 -2 0 -2 -1 0 1 2 x x q > 0 q = 0 -0.1 -0 .
It is seen that the topology of the emitted vorticity field is qualitatively the same in the three cases: a quadripolar pattern made of two counter-rotating vortex pairs is generated. This can be qualitatively interpreted as the result of a baroclinic effect of the form -(∇p × ∇ρ)/ρ 2 , in which the positive pressure gradient is related to the pressure jump across the shock wave. From that expression, it is seen that the case of a light disturbance with a negative amplitude parameter ǫ would lead to a vorticity pattern with opposite sign, i.e. a pattern made of four vortices rotating in the opposite sense to those found for a heavy density spot.

The main effects of the heat source term being i) an amplification (resp. damping) of the amplitude of the emitted perturbations and ii) an increase (resp. decrease) of the anisotropy of the emitted pattern for endothermic (resp. exothermic) case, when compared to the adiabatic case. In the strong endothermic case the amplitude of the four vortices are nearly equal, while the second vortex pair is weaker in other cases. This is consistent with the fact that the effective shock-induced compressive effect is stronger in the endothermic case, as observed in Section 2.

The emitted acoustic field is illustrated here in the bottom plots of Fig. 8, in which the acoustic pressure field p ′ γ p is plotted at time t = 4a c2 . A more global view at the interaction physics is obtained looking at the energy of the emitted waves along with the part associated to each Kovasznay mode, according to the extended definition derived in Section 4. The area used to compute the sum in Eq.(4.13) is taken equal to 12D × 12D, which was checked to be large enough to get fully converged values, with D defined as the radius of the incident Gaussian spot (see Appendix).

Results in the (M 1 , q) plane normalized by the energy of the incident density spot are displayed in Fig. 9 for the far-field solution, i.e. the transient contribution of acoustic non-propagative waves is omitted. Profiles along the q = 0 and the M 1 = 2 lines are also shown in Fig. 10.

It is observed that the total emitted energy is an increasing function of the incoming Mach number M 1 , and that the respective importance of each mode is strongly influenced by the heat source term q. In the neutral case q = 0, the emitted energy is mainly due E y and E v , i.e. to the concentration mode and the vorticity mode, the former being dominant for M 1 < 4. It is worth noting that the energy of all emitted modes is an increasing function of M 1 , excepted E s which decreases for 1 M 1 2.6 Varying q at fixed M 1 makes a more complex behavior to appear. The emitted energy is mostly related to the vorticity mode in the endothermic case, the solution being dominated by the concentration mode for sligthly negative q and exothermic cases. This is due to the case that the concentration mode is the only one which exhibit an increase for increasing Figure 10: Energy of the emitted disturbances in the case of an incident Gaussian density spot versus M 1 for the adiabatic case (q = 0) and versus q for M 1 = 2. Total energy and the part associated to each Kovasznay mode are displayed. E tot (solid thick line), E y (solid), E v (dotted), E a (dashed), E s (dotted-dashed).

q, while a decrease of the total emitted energy E tot associated to a monotonic decrease of all other modes is observed. A very fast decrease of E s is observed, leading to the fact that the entropy mode is very strong in the highly endothermic case, while it is the weakest mode in the neutral and and exothermic cases.

Gaussian entropy spot

This Section is dedicated to the interaction with a Gaussian entropy spot, and therefore is an extension of the previous analysis provided in [START_REF] Fabre | Linear interaction of a cylindrical entropy spot with a shock[END_REF] for the adiabatic case q = 0. The upstream entropy spot is defined by setting G ′ = s ′ cp in (6.2). The emitted entropy far field, vorticity far field and acoustic pressure far field are displayed in Fig. 11. The emitted disturbance topology is the same as in the density case: a quadrupolar pattern made of two counter-rotating vortex pairs is generated downstream the shock, whose intensity and anisotropy are decreasing functions of the heat source term q. The key mechanisms for vorticity generation can again be interpreted as a kind of baroclinic production term associated to the pressure jump across the shock and the density gradient associated to the entropy disturbance, see Eq. (3.4).

The total energy of the emitted far-field solution (normalized by the energy of the incident spot) and the part associated to each Kovasznay component are plotted in Fig. 12 in the (M 1 , q) plane, while profiles along the M 1 = 2 and q = 0 lines are shown in Fig. 13. It is worth noting that the concentration mode energy remains null downstream the shock, i.e. E y = 0, since it is null upstream the shock and the the concentration fluctuation is continuous at the shock according to Eq. (5.1). Some interesting differences with the density spot case are observed, which are due to the fact that the entropy spot combines a density disturbance and a temperature disturbance. First, in the adiabatic case q = 0, the normalized total emitted energy is not a monotonous function of the upstream Mach number M 1 . A decrease is observed for M 1 < M crit ≃ 2.7, which is due to a decrease of the energy of the emitted entropy mode, which is a monotonic decaying function of M 1 . The emitted acoustic and vorticity energy component, E a and E v , are growing with M 1 , E a being negligible in all cases. Therefore, the emitted field is dominated by the entropy mode for M 1 < M crit , while the vorticity mode is dominant at higher Mach number. This picture is very different from the one observed for the density spot, and it it stable with respect to a change in the parameter q. Here, the energy of all emitted modes decays when increasing q, including the emitted vortical energy which was an increasing function of q in the density spot case.

Gaussian vorticity spot

The last case deals with the interaction between a planar shock wave and a Gaussian vorticity spot, which is a model of a weak vortex. The shock/vortex interaction has been addressed by several authors, mainly via Direct Numerical Simulation, but the present analysis is the first one to cover the full (M 1 , q) plane within the LIA framework.

Results for the emitted non-acoustic fields are shown in the first 4 plots of Fig. 14. The concentration field remains uniform, as in the case of the entropy spot. A first observation is that the topology of the emitted field is different from the one observed for both incident density and entropy spot. As a matter of fact, while two vortex pairs with variable intensity were found previously, the present field is made of a strong counterrotating vortex pair, with two companion pairs of much weaker vortical structures.

The topology of the downstream acoustic field is investigated in the bottom plots of Fig. 14 which displays the generated pressure. A compression wave followed by a dilatation wave is observed, while in the two other cases the dilatation wave is emitted first.

The energy of the emitted field split into model components, normalized by the energy of the incident spot, is displayed in Figs. 15 and 16. It is observed that, in all cases, the emitted energy is dominated by the vortical component. In the adiabatic case, the acoustic energy remains larger than the entropy mode energy at all Mach number. The Figure 13: Energy of the emitted disturbances in the case of an incident Gaussian entropy spot versus M 1 for in the adiabatic case (q = 0) (top) and versus q for M 1 = 2 (bottom). Total energy and the part associated to each Kovasznay mode are displayed. E tot (solid thick line), E v (dotted), E a (dashed), E s (dotted-dashed).

opposite trend can be observed in strongly endothermic cases. All energy components are growing functions of M 1 and decreasing functions of q.

Optimal mixed disturbances with minimal radiated noise

The purpose of this section is to illustrate the possibility of finding upstream disturbances associated with peculiar emitted field. To this end, it is chosen to find the optimal combination of the three above elementary spots for minimal radiated noise.

Let us identify the emitted pressure perturbation as p ′ Y , p ′ s and p ′ ω for the density, entropy and vorticity gaussian elementary spots. Next, we introduce Π, the radiated noise emitted through the shock as

x x q > 0 q = 0 -25 -7 -7 -7 1 1 1 1 1 1 1 4 4 4 4 0 2 4 -2 -1 0 1 2 -6 -3.5 -1 .5 -1 .5 0.2 0.2 0.2 0.2 0.2 0.2 0. 6 0. 6 0.6 1 -4 -2 0 -2 -1 0 1 2 x x q < 0 aω ′ ū -4 -3 -1 -1 0.2 0.2 0. 2 0.2 0.2 0.4 0.4 0. 6 0 2 4 -2 -1 0 1 2 -6 -3.5 -1 .5 -1 .5 0.2 0.2 0.2 0.2 0.2 0.2 0. 6 0. 6 0.6 1 -4 -2 0 -2 -1 0 1 2 x x q > 0 q = 0 -0 .0 2 -0 .
Π(a s , a Y ) = V p ′ ω + a s p ′ s + a Y p ′ Y γ p 2 dx dy, (6.3)
where V is the volume of fluid after the shock. This corresponds to the acoustic perturbation obtained through combination of elementary spots as illustrated in Fig. 7, with coefficients (a Y , a s , 1) for the three elementary spots. These combination can be interpreted as a family of low-density hot vortices.

Figure 17 presents the result of the minimization of Π (6.3). The top two plots show the normalized relative amplitudes 6.4) having found that a Y < 0 over the explored range of (M 1 , q). Note that the opposite sign found for a s and a y found to minimize Π could have been intuited from Figs. 8 and 11, the density and entropy spots leading to relatively similar emission patterns. The bottom two plots compare Π for the optimal (a s , a y ) with Π(0, 0), the noise radiated by the elementary gaussian vorticity spot, showing that the vortex emitted noise was reduced by 80 to 90% by superimposing the adequate density and entropy perturbations. Figure 18's top plot shows the resulting pressure field in the case M 1 = 2 and q = 0, for which we found a Y = -0.976 and a s = 2.407. It is obtained through linear combination of the emitted pressure for the elementary spots of Figs 8, 11 and 14 with weights (a Y , a s , 1). From the levels of the emitted pressure, it is clear that the radiated noise is significantly reduced compared to either elementary spot -by 82.6%, as seen in Fig. 17. The bottom plot of Fig. 18 shows the vorticity pattern downstream of the shock for the same perturbation, following Fig. 7. Following the above procedure, it is straight-forward to minimize other fluctuations, such as vorticity, temperature, etc. Figure 17: Relative spots amplitudes minimizing the radiated noise Π (top) dependence with M 1 for q = 0 (left), and with q for M 1 = 2 (right). The corresponding integral radiated noise Π (6.3) is shown for the elementary vortex Π(0, 0) (dashed) and the optimal combination Π(a s , a Y ) (solid) in the bottom two plots. The dot-dashed line, plotted in the right axes, shows the noise reduction Π(as,aY )-Π(0,0) Π(0,0)

α Y = -a Y 1 -a Y + a s , α s = a s 1 -a Y + a s , α ω = 1 1 -a Y + a s , (
.

Concluding remarks

A complete LIA framework for the interaction between a planar shock and a Gaussian disturbance including thermal effects at the shock front was proposed, along with adequate extension of the energy of the disturbances. General expressions for the emitted field are also provided, allowing for a straightforward reconstruction of the solution. Such a framework can provide a deep insight into shock/mixed disturbances interaction, but also very acurate benchmark solutions for numerical scheme validation. Another results is the extension of Chu's definition of disturbance energy to the present framework, leading to a mathematically-grounded meaningful definition of the energy of both upstream and downstream fields. It is worth noting that mixed solutions based on the combination of the three elementary solutions analyzed in the previous section can also be very easily obtained by linear combinations of the instantaneous elementary fields. This way, some solutions with peculiar features can be obtained. This is illustrated by the search of upstream vortex-like disturbances with minimal emitted pressure perturbations. In a similar way, combining an heavy density spot with a cold entropy spot one can obtain an emitted field with a very small residual vorticity. Solutions that minimize or maximize the energy of a given emitted Kovasznay mode can be obtained, the relative weight of each upstream mode being a function of the upstream Mach number and the heat source parameter q.

Strong of a wide variety of covered shock/spot interaction configurations, this work may serve as benchmark for the development of shock-capturing numerical methods.

Appendix: Mathematical formulation for Gaussian spot/shock interaction

One addresses here the general formulation for the upstream and downstream field associated to the interaction between a 2D Gaussian non-acoustic disturbance and a planar shock within the Kovasznay decomposition framework. The upstream disturbances are advected at the constant speed U , and are then frozen in the reference frame moving with the base flow far from the shock wave.

Considering an upstream perturbation of Kovasznay mode i of the form (in cylindrical coordinates (r, φ)) f i (r) = ǫe -r 2 /D 2 , (8.1)

where D denotes the radius (taken equal to 1) of the Gaussian profile and f i the fluctuation of the physical quantity associated to the Kovasznay mode (e.g. density or specific volume for a density spot, entropy for an entropy spot, vorticity for a vortex), its decomposition into Fourier modes obtained via a 2D polar Fourier transform is (see [START_REF] Fabre | Linear interaction of a cylindrical entropy spot with a shock[END_REF] for a full description of intermediary algebra, with (k, α) the polar coordinates in the Fourier space) (8.4)

The complete emitted field is given by

s ′ 2 C p2 = ǫ π π 2 α=-π 2 Re(Z is )K r s c -Im(Z is )L r s c dα, (8.5) ω ′ 2 U 2 = - ǫ π π 2 α=-π 2 1 c Re(Z iv )K ′ r s c -Im(Z iv )L ′ r s c dα, (8.6) p ′ 2 γP 2 = ǫ 2π Re π 2 α=-π 2 ∞ k=0
Z ia e -k 2 4

e -kaη(x+M2c2t) e ika.x k dk dα , (8.7) e -kaη(x+M2c2t) e ika.x k dk dα , (8.9)

T ′ 2 T 2 = ǫ π π 2 α=-π 2 Re(Z is )K r s c -Im(Z is )L r s c dα + (γ -1) p ′
u ′ r2 U 2 = - ǫ π π 2 α=-π 2 Re(Z iv )K r s c -Im(Z iv )L r s c cos α s dα + ǫ 2π Re π 2 α=-π 2 ∞ k=0 Z ia sin α a M 2 ζ e -k 2 4
e -kaη(x+M2c2t) e ika.x k dk dα , (8.10) (8.11) where it is recalled that k a .x = k a cos α a (x + U 2 t) + sin α a (y -U 2 t cot α s ) . The auxiliary functions L, K ′ , L ′ , r s and c are defined as :

u ′ φ2 U 2 = ǫ π π 2 α=-π 2 Re(Z it )K r s c -Im(Z it )L r s c dα,
L(z) = ∞ k=0 k 2 e -( k 2 )
2 sin(kz)dk = √ πze -z 2 , (8.12)

K ′ (z) = -2z - √ π(1 - 2z
2 )e -z 2 erfi(z), (8.13)

Figure 1 :

 1 Figure 1: Downstream Mach number M 2 (left) and compression factor m (right) as functions of the upstream Mach number M 1 and the heat source/sink parameter q according to Eqs.(2.6) to (2.10).

Figure 3 :

 3 Figure3: Isovalues of the critical angle α c as a function of the upstream Mach number M 1 and the heat source term q. White areas correspond to unphysical configurations that violate the realizability constraint on the downstream flow.

  from Eq. (3.10) with β = α α < αc, (αa, η) from Eq. (3.13) with β = α α > αc, (αa, η) from Eq. (3.14) with β = α acoustic perturbation

Figure 4 :

 4 Figure4: Emission angles and attenuation factor η, obtained following the procedure summarized in Tab. 2, for γ = 1.4, M 1 = 2 and q = -2.25. Plain line: α a , dashed-line: α s , dotted-line: attenuation η. In the left plot, the additional dot-dashed line represents the non-linear dependence of β as a function of α (3.7) in the case of an acoustic incident perturbation.

Figure 5 :

 5 Figure 5: Real part (plain line) and imaginary part (dashed) of Z ii as a function of the incident wave angle α, for γ = 1.4, M 1 = 2 and q = -2.25. The corresponding emitted wave vectors angles α a and α s are those represented in Fig 4.

Figure 6 :

 6 Figure 6: Real part (plain line) and imaginary part (dashed) of Z Y i as a function of the incident wave angle α, for different incident mass fraction wave: (A r t , A Cv t ) = (1,0) (top), (0,1) (bottom). The remaining parameters are identical to Fig. 5: γ = 1.4, M 1 = 2 and q = -2.25. The corresponding emitted wave vector angle α s can be found in Fig 4.

Figure 7 :

 7 Figure 7: Emitted vorticity for incident gaussian density Y , entropy s, and vorticity ω spots. The fourth spot corresponds to the sum Σ of the three gaussian spots, resulting in yet another vorticity pattern. The dashed line illustrates the corrugated shock.

Figure 8 :

 8 Figure 8: Incident Gaussian density spot: emitted entropy, vorticity and acoustic perturbations (from top to bottom). Left: adiabatic vs endothermic case. Right: adiabatic vs exothermic.

Figure 9 :

 9 Figure 9: Energy of the emitted disturbances in the case of an incident Gaussian density spot in the (M 1 , q) plane. Total energy E tot and the part associated to each Kovasznay mode are displayed, with E y : energy of the concentration mode; E v : energy of the vorticity mode; E a : energy of the acoustic mode; E s : energy of the entropy mode.

Figure 11 :

 11 Figure 11: Incident Gaussian entropy spot: emitted entropy, vorticity and acoustic perturbations (from top to bottom). Left: adiabatic vs endothermic case. Right: adiabatic vs exothermic.

Figure 12 :

 12 Figure 12: Energy of the emitted disturbances in the case of an incident Gaussian entropy spot in the (M 1 , q) plane. Total energy E tot and the part associated to each Kovasznay mode are displayed, with E v : energy of the vorticity mode; E a : energy of the acoustic mode; E s : energy of the entropy mode.

Figure 14 :

 14 Figure 14: Incident weak vortex/Gaussian vorticity spot: emitted entropy, vorticity and acoustic perturbations (from top to bottom). Left: adiabatic vs endothermic case. Right: adiabatic vs exothermic.

Figure 15 :Figure 16 :

 1516 Figure 15: Energy of the emitted disturbances in the case of an incident weak vortex/Gaussian vorticity spot in the (M 1 , q) plane. Total energy E tot and the part associated to each Kovasznay mode are displayed, with E v : energy of the vorticity mode; E a : energy of the acoustic mode; E s : energy of the entropy mode.

Figure 18 :

 18 Figure18: Emitted pressure (top) and vorticity (bottom) fields for the optimized vortex perturbation minimizing the radiated noise Π, in the case M 1 = 2 and q = 0.

  cos(α -φ).

Table 1 :

 1 Summary of the LIA literature. ω, s, p and Y indicate the considered incident Kovásznay modes, and ∆Q the presence of a heat releasing and/or absorbing shock. T urb

† Email address for correspondence: pierre.boivin@univ-amu.fr

Acknowledgements

This work was supported by Labex MEC (ANR-10-LABX-0092) and the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the "Investissements d'Avenir".

The emitted density field ρ ′ 2 and specific volume τ ′ 2 = 1/ρ ′ 2 are computed for each component from T ′ , Y ′ and p ′ in a straightforward way using the linearized equation of state, i.e. (8.16)


Michael Baumann
email: michael.baumann@imag.fr

Vincent Daanen

Antoine Leroy

Jocelyne Troccaz

3-D Ultrasound Probe Calibration for Computer-Guided Diagnosis and Therapy

published or not. The documents may come

Introduction

Until recently, 3-D ultrasound (US) volumes had to be manually reconstructed from a number of 2-D US slices acquired while slowly moving a 2-D probe over the target region. The so-called 3D freehand method is time-consuming, imprecise and not usable for many clinical applications requiring real-time acquisition [START_REF] Rohling | A comparison of freehand three-dimensional ultrasound reconstruction techniques[END_REF]. The emergence of 3-D swept-volume US probes solved most of the enumerated problems: a mechanical device capable of sweeping the 2-D crystal array of the probe over a target region makes it possible to acquire 3-D US volumes accurately and almost in real-time (1s to 4s per acquisition) 3 .

These new capabilities open an entire new field of applications in the domain of computer guided medical interventions based on US imaging. One can imagine tool guidance systems that would operate with permanently updated US volumes, visualizing for instance slices at the tool tip position. More sophisticated applications could carry out target localization inside the volumes through real-time registration and segmentation techniques, thus allowing to match pre-operative planning with intra-operative data.

However, US-based guidance often requires knowledge about the position and orientation of the US volume in space. When using a tracking system this can be achieved by calibrating the US acquisition volume with a tracking reference fixed on the probe. Unfortunately it is virtually impossible to derive the calibration parameters directly from the geometry and parameterization of the probe. Almost all existing calibration systems rely therefore on statistical or segmentation-based object matching methods.

Calibration Methods Overview

A variety of techniques for 2-D US calibration was proposed in the literature; a comprehensive review being given in [START_REF] Mercier | A review of calibration techniques for freehand 3-D ultrasound systems[END_REF]. Calibration methods can be classified with respect to the target geometry they rely on. Single-point target methods identify a point, i.e. a bead, a calibrated pointer tip or a cross-wire, in the US image [START_REF] State | Case study: Observing a volume rendered fetus within a pregnant patient[END_REF][START_REF] Prager | Rapid calibration for 3-D freehand ultrasound[END_REF][START_REF] Muratore | Beam calibration without a phantom for creating a 3-D freehand ultrasound system[END_REF]. The difficulties consist in automatic geometry extraction in the US slice and US beam alignment with the phantom.

Multi-point target phantoms are extensions of the single-point bead or cross-wire phantoms. They consist of a number of point targets with precisely known coordinates in phantom space. Their geometric configuration makes it possible to derive the calibration parameters from the distances between the reconstructed intersection points visible in the 2-D US scan [START_REF] Leotta | An efficient calibration method for freehand 3-D ultrasound imaging systems[END_REF]. Compared to single-point phantoms they require less image acquisitions due to their more discriminative geometry, but share the phantom alignment and feature extraction problems.

Z-fiducial or N-fiducial phantoms address the alignment problem of point target methods. A calibration point is determined from the intersection points of a number of nylon strings with the US beam. This is possible due to a sufficiently discriminative wire geometry [START_REF] Comeau | Integrated mr and ultrasound imaging for improved image guidance in neurosurgery[END_REF][START_REF] Pagoulatos | A fast calibration method for 3-D tracking of ultrasound images using a spatial localizer[END_REF][START_REF] Linseth | Probe calibration for freehand 3-D ultrasound[END_REF]. Fiducial methods are more robust than point target methods but the difficulties concerning fully automatic feature extraction subsist. Also, Z-or N-fiducial phantoms require a high manufacturing accuracy to achieve a satisfying calibration quality.

Wall phantom methods are based on detection of the intersecting line of a planar surface with the 2-D US beam. In [START_REF] Prager | Rapid calibration for 3-D freehand ultrasound[END_REF], a water tank bottom is imaged for calibration. The authors of [START_REF]Ultrasound guided surgery: Image processing and navigation[END_REF] address the reverberation and line thickness problems inherent of wall phantoms by using a membrane variant. Both phantoms have difficulties when confronted with steep angles between the US beam and the plane because they cause line intensity and line sharpness degradation [START_REF] Mercier | A review of calibration techniques for freehand 3-D ultrasound systems[END_REF]. The Cambridge phantom scans a rotating bar, thus creating a virtual plane, to solve these problems [START_REF] Prager | Rapid calibration for 3-D freehand ultrasound[END_REF]. The advantage of plane phantoms lies in the robustness of the feature extraction process which can, as a consequence, be reliably and fully automated. The pitfall of this method lies in the nondiscriminative phantom geometry which can result in underdetermined systems if the acquired calibration samples do not cover all degrees of freedom. This can be avoided by strictly respecting the acquisition protocols presented in [START_REF] Prager | Rapid calibration for 3-D freehand ultrasound[END_REF][START_REF] Treece | High-definition freehand 3-D ultrasound[END_REF].

Registration Phantoms: the last class of calibration methods relies on surface or intensity based registration techniques and therefore has the advantage of being inde-pendent of phantom geometry. The only requirement on phantom shape is that its US image is sufficiently discriminative with respect to rotations and translations, which is true for non-symmetric phantoms. The lack of precision of registration algorithms is the major drawback of this approach. To our knowledge, only one study examined registration-based 2-D probe calibration, registering US slices with an MRI image of the phantom [START_REF] Blackall | An image registration approach to automated calibration for freehand 3-D ultrasound[END_REF]. A 3-D approach is discussed in the next paragraph.

3-D Probe Calibration

Until today, only few studies about calibration of 3-D probes were carried out. Poon and Rohling [START_REF] Poon | Comparison of calibration methods for spatial tracking of a 3-D ultrasound probe[END_REF] compared 3-D calibration based on a IXI-fiducial wire phantom, a pointer tip phantom and a cube phantom. The IXI wire phantom and the cube phantom methods require only one volume acquisition for calibration. The presented feature detection is semi-automatic. The best results yielded the IXI phantom with a mean error in reproducibility of 1.5 mm, a RMS error of the point accuracy measure of 2.15 mm and a RMS error of the reconstruction accuracy by distance measure of 1.52 mm. Bouchet et al [START_REF] Bouchet | Calibration of three-dimensional ultrasound images for image-guided radiation therapy[END_REF] examined Z-fiducial phantom and achieved a RMS point accuracy error of 1.1mm. Two variants of a surface registration based 3-D calibration method were presented by Lange and Eulenstein in [START_REF] Lange | Calibration of swept-volume 3-D ultrasound[END_REF]. The first one registers 3-D US images of the phantom with a geometric model derived from its CT scan. The second variant registers a number of US images of the phantom acquired from different positions. In both cases, surfaces are extracted manually. The authors claim that the latter approach could be fully automated. The CT variant performed best and yielded a RMS error in reproducibility precision of 1.8 mm and a RMS error in point accuracy of 2.0 mm. The ultrasound speed distortion problem is not addressed.

In this study we propose a 3-D US calibration method based on a single plane membrane phantom. A fast, precise and accurate 3-D feature extraction algorithm relying on the 2-D Hough transform is presented. In contrast to existing 3-D US calibration systems, the feature extraction process is fully automated. In the result section, precision and accuracy assessments are carried out using a specially designed validation phantom.

Materials and Methods

Acquisition Hardware

The acquisition hardware consists of a GE Voluson 730 Pro 3-D US scanner and a NDI Polaris optical tracker with a 0.25 mm RMS error (as communicated by NDI). The tracking system operates with wireless (passive) infrared-reflecting rigid bodies equipped with flat markers. The ultrasound volumes are acquired with a 5 to 9 Mhz two dimensional curved array probe (see Fig. 1a). The piezo array of the probe is mounted on a mechanical device which is capable of sweeping regularly around its rotation axis within a predefined angular range. During the continuous sweeping process the US hardware reconstructs 3-D volumes from the series of acquired 2-D slices. The 3-D acquisition time ranges from 1s to 4s, depending mainly on sweep angle and axial acquisition depth. Images are digitally transferred using a proprietary software from GE Medical Systems named 4D View. The US scanner also communicates the voxel size. The scan converter assumes the speed of sound (SoS) in tissue to be 1540 m/s.

The Membrane Phantom

The calibration phantom being dedicated to a clinical context, ergonomics considerations had an important impact on its design. We use a variation of the wall phantom presented in [START_REF] Prager | Rapid calibration for 3-D freehand ultrasound[END_REF], which is based on imaging the bottom wall of a water tank. The geometric form of the wall, which is a line in 2-D and a plane in 3-D, can be very robustly extracted from the US data using statistical algorithms like the Hough transform. This makes it possible to fully automate the feature extraction process without significant loss in precision and accuracy. This represents a big advantage over semi-automatic point-detection based phantoms in terms of calibration speed and ease of use. To overcome the plane thickness and the reverberation problems observable in US images of rigid surfaces [START_REF] Prager | Rapid calibration for 3-D freehand ultrasound[END_REF] a filigrane nylon mesh membrane, tightly spanned on a planar rigid support with a circular and about 20cm wide hole, is used as target (see Fig. 1b). Reverberation is further reduced by inclining the membrane plane with respect to the water tank bottom by 45 • . A tracking reference (rigid body) is mounted on the membrane frame for phantom localization. The phantom is filled with water and equipped with a thermometer to measure water temperature.

3-D Calibration Mathematics

As illustrated in Fig. 2, four references are relevant for calibration: first of all, the membrane space M is defined as a reference in which the membrane lies in the origin and is parallel to e x and e y base vectors of M. In this space, every point with a zero z-ordinate is a membrane plane point. The phantom space Ph and the probe space Pr are defined by the rigid bodies that are attached on the phantom and on the probe. Finally, the US volume space U corresponds to the voxel space of the 3-D images acquired by the ultrasound device. T Ph2M , T Pr2Ph and T U2Pr are homogenous 4x4 transformation matrices.

Suppose that we identified a point p = (x, y, z) in a US volume U as a point belonging to the membrane. With s = (s x , s y , s z ) denoting the voxel scale factors, it is verified that

    m 1 m 2 0 1     = T Ph2M • T Pr2Ph • T U2Pr •     s x x s y y s z z 1     .
(

) 1
where T Ph2M is known from membrane pre-calibration (see chap. 2.4) and T Pr2Ph is given by the tracking system. Further, the scaling vector s is communicated by the US hardware. The remaining unknown element is the homogenous rigid transformation T U2P . For convenience we define the elements of T Pr2Ph (

) 2
Using Euler angles and a three-dimensional vector we can represent T U2P with six variables, which leaves us in total with 6 unknowns to solve for. A detected plane can be added to the equation system by adding at least three plane points (Using of course the T Pr2Ph measured while acquiring the US volume in which the plane was detected).

Fig. 2. Illustration of the transformations involved in the calibration process. Note that scaling is omitted from the scheme for simplification.

Membrane Pre-calibration

To reduce the number of degrees of freedom of the calibration process, the membrane space M is determined using a pointer equipped with a tracking reference. A large number of surface points of the membrane-supporting structure is acquired in order to compute the plane equation using a least square approximation combined with a simple M-estimator to increase robustness. Since the phantom rigid body is permanently fixed on the phantom, pre-calibration has to be carried out only once.

Acquisition Protocol

The major drawback of single-plane phantoms resides in their barely discriminative geometry. A plane can be described with only three variables, from which it follows that even with an optimal acquisition protocol, a minimum number of two acquisitions is necessary to cover all degrees of freedom. To obtain robust results the twelve-step acquisition protocol presented in [START_REF] Prager | Rapid calibration for 3-D freehand ultrasound[END_REF] is used. The protocol improvement presented in [START_REF] Treece | High-definition freehand 3-D ultrasound[END_REF] mainly addresses the z-axis imprecision problem inherent of most 2-D calibration systems. Since 3-D probes give as much information in z-direction as in x-or y-direction this modification yields no particular advantage in the 3-D domain but requires at least 18 steps. For that reason we stick to the original version. Sweeping and volume reconstruction being a continuous process of 1 to 4 seconds, significant distortions can be observed in the US volume when the probe is moved rapidly. Also, no direct access to the digital data is available which prevented synchronization of probe position measurement with US image acquisition. Therefore an articulated arm for complete probe immobilization during acquisition is used, eliminating all motion-induced artifacts and time lags. Furthermore, immobilizing the probe makes it possible to perform high precision position measuring based on a large number of measures and outlier elimination.

Feature Extraction

The first step of the feature extraction process consists in correcting the distortion caused by the difference between US speed in water at room temperature and in human tissue at 37 • . To determine US speed in water in function of temperature the polynomial formula established by Bilaniuk and Wong was used [START_REF] Bilaniuk | Speed of sound in pure water as a function of temperature[END_REF][START_REF] Bilaniuk | Erratum: Speed of sound in pure water as a function of temperature[END_REF]. A distortion geometry overview for all common probe types is given in [START_REF] Goldstein | The effect of acoustic velocity on phantom measurements[END_REF]. The distortion geometry of a sectorial probe is given in Fig. 3a. With v t W being the US velocity in water for a given temperature t and v T being the velocity in tissue, d T is determined using the following formula:

d T = v T v W • d W . (3)
Sectorial probe speed correction requires manual definition of the US origin and the scan head surface radius. A graphical user interface was developed for this purpose (see Fig. 3b). Origin and surface radius have to be defined only once during calibration. Plane detection can be carried out using the 3-D Hough transformation, but it would take several minutes to compute the result. Fortunately it is possible to determine the plane with good precision by simply extracting its intersection with two arbitrary volume slices, using the 2-D Hough transform. To facilitate and to accelerate US speed correction the xy and zy planes passing through the scan head origin were used. The Hough transform implementation uses intensity accumulation and the following threshold s H for an image I:

s H = max{i ∈ Hist(I)} + (max{i ∈ I} -min{i ∈ I})/3 . ( 4
)
The purpose of s H is to ignore the low-intensity water background, which represents the largest part of the image.

Optimization

Optimization of ( 2) is carried out with the non-linear Levenberg-Marquardt implementation given in [START_REF] Press | Numerical Recipies in C, 2nd edition[END_REF]. A random restart scheme within a range of reasonable initialization values robustifies this process.

Visual Back-tests

The plane coordinates as resulting from the optimization process are visualized as a line in the slices used for feature extraction, and distance plus rotational errors between the segmented and the calculated line are evaluated. This allows for manual replacement of evident outliers with new acquisitions and for recomputation of the calibration without requiring a complete restart. 3 Experimental Results

Test Configuration

Precision and accuracy assessments were carried out using a membrane plane precalibration with an RMS surface distance error of 0.43mm for the measured surface points. A total of ten calibrations were performed using the twelve-step protocol. The probe rigid body was not moved between calibrations. The probe was mounted on an articulated arm and immobilized during position and image acquisition. The water temperature was 23 • . The acquired US volumes had a size of (199, 199, 199) isotropic voxels with 0.477mm side lengths.

Feature Extraction Quality

The Hough transform extracted lines correctly for 238 out of 240 acquired images. Line detection failed for images on which only a very small part of the membrane was visible. In these cases, lines had to be manually determined. Note that no manual outlier elimination was carried out. As detection failures were rare, calibration was in average carried out in about twenty minutes. To get a better idea about the quality of feature extraction and about the presence of distortions in the membrane images we measured the detection precision: using (2) we can calculate the distance error between a measured plane point and the computed plane as follows:

ǫ(x,
For each calibration, the average and the root mean square (RMS) distance of a set of points to the pre-calibration plane was computed using [START_REF] Muratore | Beam calibration without a phantom for creating a 3-D freehand ultrasound system[END_REF]. For each line we computed ten equidistant points between the extreme points on the line segment inside the US volume. The angular feature extraction error is defined as the angle between the computed plane normal and the cross product of the directional vectors of the two extracted lines.

Based on this definition the maximum and the RMS angular errors were computed for each acquired volume of the calibration. The aggregated errors for all calibrations can be found in Table 1.

Calibration Precision

The calibration precision measures the reproducibility of calibration results. Again, both the translational and the angular errors were assessed. The translational error is defined as the standard deviation of the volume center after scaling and right-hand multiplication to the different calibration transformations T i U2Pr . The angular error is measured as the standard deviation of angular differences between the (0, 0, 1) vector after scaling and right-hand multiplication to the different calibration transformations T i U2Pr (see Table 2).

Reconstruction Accuracy

Reconstruction accuracy was assessed using the bead phantom illustrated in Fig. 5. Note that the beads are co-planar within a precision of 0.25mm (RMS). The left-hand three beads form the left triangle while the right-hand beads form the right triangle.

The distance d B of the triangle barycenters was evaluated with an estimated accuracy of about 0.5mm. Twenty US images of the phantom were acquired, ten imaging the left triangle and ten the right triangle. Note that images of one triangle did not intersect with the second one. The bead centers were manually extracted from the images. The reconstructed triangle barycenters and normals were then projected into probe space for each calibration, which yielded 100 point and vector pairs. The distance error for a point pair is defined as the difference between their Euclidean distance and d B . The angular error for a vector pair is defined as the angle between both vectors. The results are given in Table 3.

Discussion

Probe calibration for currently available swept-volume 3-D probes makes only sense for static applications. Real-time access to acquired 3D images is currently not provided by the Voluson hardware. Also, depending on scan parameters, the duration of the volume sweeping process ranges from 1s to 4s. Due to probe or tissue motion, the physical location of voxels can therefore be way off the position indicated by the calibration. The latter problem could be reduced by scan-heads equipped with high-frequency sweepingdevices, but it will disappear only with non-sweeping 2-D piezo array US probes. Until then, it is better to calibrate probes using an articulated arm for immobilization. Passing from 2-D to 3-D calibration improves calibration results (for information on 2-D precision of the membrane phantom see [START_REF] Hook | Probe calibration for 3-D image localization[END_REF]) because the z-axis-uncertainty inherent of 2-D calibration is eliminated: when giving a plane instead of a line as input to the optimizer the rotational degrees of freedom are significantly better covered. This allowed us to reduce the number of acquisitions for calibration while still achieving precise and accurate results.

Feature extraction from membrane phantom images showed both robust and precise results. On our set of test images the line blurring and intensity degradation effects occurring when scanning a wall phantom from an oblique angle were correctly handled by the Hough transform: lines were consistently placed in the center of the beam width. Note that the feature extraction precision RMS and maximum errors reported in Table 1 are relatively small, which indicates that the physical plane location corresponds indeed with the beam width center line. Membrane reverberation was not observable and did therefore not disturb the detection process. Also, the membrane phantom was not exposed to line thickness problems. Due to these characteristics, feature extraction could be fully automated (up to the manual US origin and probe radius determination required for US speed correction).

The user-independency resulting from automated feature extraction is partially counter-balanced by the necessity to follow a protocol for data acquisition, which can reintroduce user bias. Nevertheless we believe that it is more convenient to follow a simple acquisition protocol instead of extracting features semi-automatically from US volumes. Further we preferred to correct for US speed errors instead of requiring 50 • water or phantom fill materials that have the same US characteristics than human tissue, which makes it possible to use the phantom in a clinical context.

The overall calibration time of about 20 minutes is mostly due to manipulation of the articulated arm, high precision probe position measurement and data transfer, which requires several manual interventions. Feature extraction and optimization is computed in several seconds. In cases where the feature extraction precision evaluation shows poor results, a visual verification and eventually a correction have to be carried out, which requires some additional minutes. Significant speed-up could be achieved by automating the communication between the US scanner and the calibration computer. The presented calibration system assumes that the SoS in human tissue is uniform and that it corresponds to the SoS internally used by the US scanner, which is in general the mean SoS in tissue of 1540m/s. However, SoS varies with different types of tissue: The SoS in fat is approximately 1450m/s, in blood 1570m/s, in the brain 1541m/s and in water 1480m/s. As the US generally crosses tissue layers of different thickness and different types on its way through the body, and as the target tissue is often viewed from various positions, the in vivo accuracy of the calibration may show fluctuations of more than 5 per cents in extreme cases during an examination. Also, for some applications it would be appropriate to use a different mean SoS than the 1540m/s for calibration, but this is beyond the scope of this study. Future work will address the twelve-step acquisition protocol which contains a lot of redundancy. Also, our system does currently not provide a foolproof indicator for missing coverage of degrees of freedom. We therefore started experiments with an Eigenvalue system similar to the one presented in [START_REF] Hsu | Rapid, easy and reliable calibration for freehand 3-D ultrasound[END_REF].

A robust 3-D US probe calibration method designed for clinical usage was presented. Calibration can be carried out in about twenty minutes due to fully automatic 3-D Plane extraction based on robust and efficient 2-D line detection. The point reconstruction accuracy of our phantom can compete with previously presented 3-D phantoms: Lange and Eulenstein communicated RMS errors between 2.0mm and 2.2mm [START_REF] Lange | Calibration of swept-volume 3-D ultrasound[END_REF], Bouchet et al were confronted to 1.1mm RMS point accuracy [START_REF] Bouchet | Calibration of three-dimensional ultrasound images for image-guided radiation therapy[END_REF] while Pohn and Rohling published errors between 1.52mm for their IXI-wire, 1.59mm for the cube and 1.85mm for their stylus approach. With 0.90mm RMS point accuracy (see Table 3) we achieved slightly better results. Finally, the proposed method is temperature-independent and uses water as transmission matter which facilitates its usage.

Fig. 1 .

 1 Fig. 1. Calibration hardware. Figure (a) shows an endorectal US probe mounted on an articulated arm. Figure (b) shows the membrane phantom. Both the probe and the membrane are equipped with infra-red reflecting passive rigid bodies for tracking.

Fig. 3 .

 3 Fig. 3. US speed correction. Fig. (a) illustrates the correction geometry of a sectorial probe. O US is the probe origin, r the probe surface radius, d W the point's P distorted distance from the probe surface in water, d T the corrected distance and P ′ the corrected point. Fig. (b) shows the probe mask used to determine US origin and probe surface radius.

Fig. 4 .

 4 Fig. 4. Screenshot of a successful automatic plane extraction. Note that lines are correctly detected in spite of a degraded membrane image caused by a steep scan angle. The arrows point at the line intersections with the mask.

Fig. 5 .

 5 Fig. 5. Reconstruction Accuracy Measurement Phantom.

  y, z) = a 31 (s x xb 11 + s y yb 12 + s z zb 13 + b 14 ) + a 32 (s x xb 21 + s y yb 22 + s z zb 23 + b 24 ) + a 33 (s x xb 31 + s y yb 32 + s z zb 33 + b 34 ) + a 34 .

Table 1 .

 1 Aggregated Feature Extraction Precision.

		Distance Error Distance Error Angular Error
		[mm]	[vox]	[deg]
	RMS Error	0.37	0.77	0.26
	Max Error	1.30	2.73	1.09

Table 2 .

 2 Calibration Precision.

		Distance Error Distance Error Angular Error
		[mm]	[vox]	[deg]
	Standard Deviation	1.15	2.41	0.61
	Max Error	1.99	4.03	1.12

Table 3 .

 3 Reconstruction Accuracy Results.

		Distance Error Distance Error Angular Error
		[mm]	[vox]	[deg]
	RMS Error	0.90	1.90	1.79
	Max Error	2.44	5.11	3.03

However, most currently available systems don't yet provide a real-time data transfer interface for 3-D data. Nevertheless one can acquire so called "4-D" images of three orthogonal volume slices in real-time using a video-capture device. In the rest of this article we make abstraction of this restriction, hoping that it will disappear with the next generation of 3-D echographs.

Acknowledgements: This project is supported by PRAXIM/KOELIS, ANRT and PHRC 2003 "Prostate-echo" (from French Health Ministry) grants.