A Gramfort

G Peyré

M Cuturi

Fast Optimal Transport Averaging of Neuroimaging Data

. It leads to a smooth convex optimization problem and an algorithm with strong convergence guarantees. We illustrate the versatility of this tool and its empirical behavior on functional neuroimaging data, functional MRI and magnetoencephalography (MEG) source estimates, defined on voxel grids and triangulations of the folded cortical surface.

Introduction

Computing the average of some observations may seem like a trivial problem, yet it remains an active topic of research in mathematics, statistics and applications such as medical imaging. The problem of atlas computation from images [START_REF] Joshi | Unbiased diffeomorphic atlas construction for computational anatomy[END_REF], or meshes [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF], or the problem of group analysis from functional imaging data [START_REF] Thirion | Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses[END_REF] are particularly relevant for this field. The challenge is that natural phenomena are usually described in terms of physical and temporal event locations, along with their intensity. While Euclidean averaging is standard and has some benefits such as low computation time, this procedure ignores the geometry of the space the observations belong to; the image of an average brain image obtained by Euclidean averaging of individual voxels does not yield the image of the brain of an average individual.

Starting from observations defined on a regular or irregular grid, our aim is to provide a model-free approach to average them that only builds upon geometric arguments. An example of such data are functional MRI (fMRI) data defined on a voxel grid or a triangulated cortical surface. The approach aims to be intrinsically geometric in the sense that it only requires the prior knowledge of a metric between the locations on the grid. The technique aims to be versatile in the sense that it can be applied to weighted samples taking values on a discretized space with no assumptions on the regularity of the metric.

The approach we propose is inspired by optimal transport theory [START_REF] Villani | Optimal transport: old and new[END_REF] and can be seen as an extension of the Wasserstein barycenter problem [START_REF] Agueh | Barycenters in the Wasserstein space[END_REF][START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF][START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF], which aims at estimating a probability measure which bests approximates a family of probability measures in the Wasserstein metric sense. The challenge of using optimal transport in the setting we consider comes from the fact that Wasserstein distances (and their barycenters) are defined for probability measures only. While some data are normalized such the Orientation Diffusion Function (ODF) in diffusion MRI [START_REF] Descoteaux | Deterministic and probabilistic tractography based on complex fibre orientation distributions[END_REF], a number of medical imaging data are non-normalized. Here we bypass this limitation using a generalization of optimal transport distances proposed by [START_REF] Kantorovich | On a space of totally additive functions, vestn[END_REF]. This extension comes at a price, since it introduces an additional parameter (the cost of adding or removing mass) which can be difficult to tune. We propose a simple way to mitigate this problem by introducing a natural constraint on the overall mass of the barycenter, which we set to be equal to the average of the masses of all samples. We provide an efficient method to compute Kantorovich means by building upon the first algorithm of [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF]. We provide intuitions on the behavior of our method and demonstrate its relevance with simulations and experimental results obtained with fMRI and MEG data which are neuroimaging data defined on a voxel grid or a triangulation of the folded cortical mantle.

This paper is organized as follows. We start in Section 2 with a reminder on optimal transport and the Kantorovich metric for non-normalized measures. We introduce Kantorovich means in Section 3 and describe efficient algorithms to compute them. Section 4 contains simulations and results on publicly available fMRI data with 20 subjects and MEG data with 16 subjects.

Optimal Transport Between (Un)Normalized Measures

We introduce in this section Wasserstein distances for non-negative measures on a finite metric space (Ω, D). Simply put, we consider normalized histograms on a grid of locations, and assume the distance between any two locations is provided. Next, we extend Wasserstein distances to non-negative vectors of bounded mass.

U (a, b) def = {T ∈ R d×d + | T 1 d = a, T T 1 d = b}. (1)
Given the constraints induced by a and b, one naturally has that U (a, b) is empty when

|a| 1 =
W p p (a, b) = OT(a, b, M p ) def = min T ∈U (a,b) T, M p , (2)
parameterized by the marginals a, b and matrix M p .

Optimal Transport for Unnormalized Measures. If the total masses of a and b differ, namely |a| 1 = |b| 1 , the definition provided above is not useful because U (a, b) = ∅. Several extensions of the OT problem have been proposed in that setting; we recall them here for the sake of completeness. In the computer vision literature, [START_REF] Rubner | The earth movers distance, multi-dimensional scaling, and color-based image retrieval[END_REF] proposed to handle that case by: (i) relaxing the equality constraints of U (a, b) to inequality constraints T 1 d ≤ a, T T 1 d ≤ b in Equation (1); (ii) adding an equality constraint on the total mass of the solution 1 T d T 1 d = min(|a| 1 , |b| 1 ); (iii) dividing the minimum of T, M under constraints (i,ii) by min(|a| 1 , |b| 1 ). This modification does not, however, result in a metric. [START_REF] Pele | A linear time histogram metric for improved sift matching[END_REF] proposed later a variant of this approach called EMD-hat that incorporates constraints (i,ii) but (iii') adds to the optimal cost T , M a constant times min(|a| 1 , |b| 1 ). When that constant is large enough M , [START_REF] Pele | A linear time histogram metric for improved sift matching[END_REF] claim that EMD-hat is a metric. We also note that [START_REF] Benamou | Numerical resolution of an unbalanced mass transport problem[END_REF] proposed a quadratic penalty between the differences of masses and made use of a dynamic formulation of the transportation problem.

Kantorovich Norms for Signed Measures. We propose to build on early contributions by Kantorovich to define a generalization of optimal transport distance for unnormalized measures, making optimal transport applicable to a wider class of problems, such as the averaging of functional imaging data. [START_REF] Kantorovich | On a space of totally additive functions, vestn[END_REF] proposed such a generalization as an intermediary result of a more general definition, the Kantorovich norm for signed measures on a compact metric space, which was itself extended to separable metric spaces by [START_REF] Hanin | An extension of the kantorovich norm[END_REF]. We summarize this idea here by simplifying it to the case of interest in this paper where Ω is a finite (of size d) probability space, in which case signed measures are equivalent to vectors in R d . [START_REF] Kantorovich | On a space of totally additive functions, vestn[END_REF] propose first a norm for vectors z in the orthogonal of 1 d (vectors z such that z T 1 d = 0), by considering the 1-Wasserstein distance between the positive and negative parts of z, z K = W 1 (z + , z -). A penalty vector ∆ ∈ R d + is then introduced to define the norm x K of any vector x as the minimal value of z K + ∆ T |z -x| when z is taken in the space of all vectors z with zero sum, and |z -x| is the element-wise absolute value of the difference of vectors z and x. For this to define a true norm in R d , ∆ must be such that

∆ i ≥ max j m ij and |∆ i -∆ j | ≤ m ij .
The distance between two arbitrary non-negative vectors a, b of different mass is then defined as a -b K . As highlighted by [26, p.108], and if we write e i for the i th vector of the canonical basis of R d , this norm is the maximal norm in R d such that for any i, j ≤ d, e i -e j K = m ij , namely the maximal norm in the space of signed measures on Ω such that the norm between two Dirac measures coincides with Ω's metric between these points.

Kantorovich Distances for Unnormalized Nonnegative Measures. [START_REF] Guittet | Extended kantorovich norms: a tool for optimization[END_REF] noticed that Kantorovich's distance between unnormalized measures can be cast as a regular optimal transport problem. Indeed, one simply needs to:

(i) add a virtual point ω to the set Ω = {1, • • • , d} whose distance D(i, ω) = D(ω, i) to any element i in Ω is set to ∆ i ;
(ii) use that point ω as a buffer when comparing two measures of different mass. The appeal of Kantorovich's formulation in the context of this work is that it boils down to a classic optimal transport problem, which can be approximated efficiently using the smoothing approach of [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] as discussed in Section 3. To simplify our analysis in the next section, we only consider non-negative vectors (histograms) a ∈ R d + such that their total mass is upper bounded by a known positive constant. This assumption alleviates the definition of our distance below, since it does not require to treat separately the cases where either

|a| 1 ≥ |b| 1 or |a| 1 < |b| 1 when comparing a, b ∈ R d + .
Note also that this assumption always holds when dealing with finite collections of data. Without loss of generality, this is equivalent to considering vectors a in R d + such that |a| 1 ≤ 1 with a simple rescaling of all vectors by that constant. We define next the Kantorovich metric on S d , where

S d = {u ∈ R d + , |u| 1 ≤ 1}. Definition 1 (Kantorovich Distances on S d ). Let ∆ ∈ R d + such that ∆ i ≥ max j m ij and |∆ i -∆ j | ≤ m ij . Let p ≥ 0.
For two elements a and b of S d , we write α = 1 -|a| 1 ≥ 0 and β = 1 -|b| 1 ≥ 0. Their p-Kantorovich distance raised to the power p is

K p p∆ (a, b) = OT([ a α ] , b β , M p ), where M = M ∆ ∆ T 0 ∈ R d+1×d+1 + . (3)
The Kantorovich distance inherits all metric properties of Wasserstein distances: the mapping which to a vector a associates a vector [a; 1 -|a| 1 ] ∈ Σ d+1 can be regarded as a feature map, to which the standard Wasserstein distance using M (which is itself a metric matrix) is applied.

Kantorovich Mean of Unnormalized Measures

Consider now a collection {b 1 , • • • , b N } of N non-negative measures on (Ω, D)
with mass upper-bounded by 1, namely N vectors in S d . Let β j = 1 -b j be the deficient mass of b j . Our goal in this section is to find, given a vector of virtual costs ∆ and an exponent p, a vector a in S d which minimizes the sum of its p-Kantorovich distances K p p∆ to all the b j , a ∈ argmin

u∈S d 1 N N j=1 K p p∆ (u, b j ) = argmin u∈S d 1 N N j=1 OT( u 1-|u|1 , b j β j , M p ). (P1)
Because of the equivalence between Kantorovich distances for points in S d and Wasserstein distances in the d + 1 simplex, this problem can be naturally cast as a Wasserstein barycenter problem [START_REF] Agueh | Barycenters in the Wasserstein space[END_REF] with metric M . Problem (P1) can be cast as a linear program with N (d + 1) 2 variables. For the applications we have in mind, where d is of the order or larger than 10 4 , solving that program is not tractable. We discuss next computational approaches to solve it efficiently.

Smooth Optimal Transport. [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF] and [START_REF] Bonneel | Sliced and radon wasserstein barycenters of measures[END_REF] have proposed efficient algorithms to solve the Wasserstein barycenter problem in low dimensional Euclidean spaces. These approaches are not, however, suitable when one considers observations on the cortex, for which all pairs shortest path metrics (inferred from a graph structure connecting all voxels) are preferred over Euclidean metrics. To solve Problem (P1) we turn instead to a recent series of algorithms proposed in [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF], [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF] and [START_REF] Cuturi | A smoothed dual approach for variational wasserstein problems[END_REF] that all exploit the regularized OT approach suggested in [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. Among these recent approaches, we propose to build in this work upon the first algorithm in [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF], which can be easily modified to incorporate constraints on a. This flexibility will prove useful in the next section. The strategy of [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF] is to regularize directly the optimal transport problem by an entropic penalty, whose weight is parameterized by a parameter λ > 0,

OT λ (a, b, M p ) def = min T ∈U (a,b) T, M p - 1 λ H(T ),
where H(T ) stands for the entropy of the matrix T seen as an element of the sim-

plex of size d 2 , H(T ) def = -ij t ij log(t ij ).
As shown by [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF], the regularized transport problem OT λ admits a unique optimal solution. As such, OT λ (a, b, M p ) is a differentiable function of a whose gradient can be recovered through the solution of the corresponding smoothed dual optimal transport. Without elaborating further on this approach, we propose to simply replace all expressions that involve an optimal transport problem OT in our formulations by their smoothed counterpart OT λ .

Sensitivity of Kantorovich Means to the Parameter ∆. The magnitude of the solution a to Problem (P1) depends directly on the virtual distance ∆. Suppose, for instance, that ∆ = ε1 d with ε arbitrarily small. In that case a should converge to a unit mass on the last (virtual) bin and would therefore be equal to the null histogram 0 d on the d other bins. If, on the contrary, ∆ = γ1 d and γ is large, we obtain that K p p∆ (a, b)/γ grows as | |a| 1 -|b| 1 |. Therefore a minimum of Problem (P1) would necessarily need to have a total mass that minimizes j |a| 1 -|b j | 1 , namely a total mass equal to the median mass of all b j . This sensitivity of the solution a to the magnitude of ∆ may be difficult to control. Choosing adequate values for ∆, namely setting the distance of the virtual point to the d other points, may also be a difficult parameter choice. To address this issue we propose to simplify our framework by introducing an equality constraint on the mass of the barycenter a in our definition, and let ∆ be any non-negative vector, typically set to a large quantile of the distribution of all pairwise distances M p ij times the vector of ones 1 d . Under these assumptions, we can now propose p-Kantorovich means with a constraint on the total mass of a. Remaining parameters in our approach are therefore only p and λ. In practice we will fix p = 1, which corresponds to the Earth Mover's Distance [START_REF] Rubner | The earth movers distance, multi-dimensional scaling, and color-based image retrieval[END_REF], and use a high λ, namely a small entropic regularization of order 1/λ, which has also the merit of making Problem (P1) strongly convex. λ is set in our experiments to 100/median(M ), where median(M ) is the median of all pairwise distances {M ij } ij .

Definition 2 (p-Kantorovich Means with Constrained Mass

). Let ∆ ∈ R d + and p ≥ 0. A Kantorovich mean with a target mass ρ ≤ 1 of a set of N histograms {b 1 , • • • , b N } in S d is the unique vector a in S d such that: a ∈ argmin a∈S d |a|1=ρ 1 N j OT λ (a, b j , M p ).
We provide in Algorithm 1 an implementation of [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF]Alg.1]. Unlike their version, we only consider a fixed step-length exponentiated gradient descent, and add a mass renormalization step. We set the default mass of the barycenter to be the mean of the masses of all histograms. We use the notation • for the elementwise (Schur) product of vectors. Note that the computations of N dual optima in line 7 of Algorithm 1 below can be vectorized and computed using only matrix-matrix products. We use GPGPUs to carry out these computations.

Application to the Averaging of Neuroimaging Data

Neuromaging data are defined on a grid of voxels, eventually restricted to the brain volume, or on a triangulation of the cortical mantle obtained by segmentation of MRI data. Examples of data most commonly analyzed on a grid are fMRI data, while neural activity estimates derived from MEG/EEG data are often restricted to the cortical surface [START_REF] Dale | Dynamic statistical parametric neurotechnique mapping: combining fMRI and MEG for high-resolution imaging of cortical activity[END_REF]. Anatomical data such as cortical thickness, which is a biomarker of certain neurodegenerative pathologies, is also defined on the surface. In all cases the data are defined on a discrete set of points (voxels or vertices) which have a natural distance given by the geometry of the brain. The Algorithm 1 p-Kantorovich Barycenter with Constrained Mass Compute all dual optima α j of OT λ (a, b j , M ) using [7, Alg.3] 8:

1: Input: {b 1 , • • • , b N } ⊂ S d , metric M , quantile q, p ≥ 0,
a ← a • exp(-c 1 N N j=1 α j ); (gradient update) 9: ai ← ρai/ d l=1 a l if i ≤ d, 1 -ρ if i = d + 1.
(projection on the simplex/mass constraint) 10: end while

11: Output a 1:d ∈ S d .
data are also non-normalized as they represent physical or statistical quantities, such as thickness in millimeters or F statistics. Such data are therefore particularly well adapted to the algorithm proposed in this paper: they are defined on discrete space, are non-normalized and there exists a natural ground metric.

The difficulty when averaging neuroimaging data is the anatomo-functional variability: every brain is different. The standard approach to compensate for this variability across subjects is to smooth the data to favor the overlap of signal of interest after the individual data have been ported to a common space using anatomical landmarks (anatomical registration). Volume data are typically redefined in MNI space while surface data are transferred to an average cortical surface, using for instance the FreeSurfer software 5 . When working with fMRI the spatial variability is commonly compensated by smoothing the data with an isotropic Gaussian kernel of Full Width at Half Maximum (FWHM) between 6 and 8 mm. MEG and EEG suffer from the same spatial variability, but also from the temporal variability of neural responses which is compensated by low pass filtering the data. By employing a transportation metric informed by the geometry of the domain, this smoothing procedure as well as the setting of the kernel bandwidth are not needed.

In the following experiments we focus on spatial averaging, although extension to spatiotemporal data is straightforward provided the metric is defined along the time axis. When working with a voxel grid the distance is the Euclidian distance taking into account the voxel size in millimeters, and when working on a cortical triangulation, the distance used is the geodesic distance computed on the folded cortical mantle. We now present results of a simulation study where standard averaging with Gaussian smoothing is compared to Kantorovich means. Simulation results are followed by experimental results obtained with fMRI data from 20 subjects and MEG data on a population of 16 subjects. Simulation setup. In this experiment using on a triangulation of the cortex, we simulated signals of interest in two brain regions using the functional parcellation provided by the FreeSurfer software. We used regions Broadman area 45 (BA45) and the visual area MT. We simulated for a group of 100 subjects random positive signals in these two regions. For each subject and each region, the signal is focal at a random location with a random amplitude generated with a truncated Gaussian distribution (mean 5, std. dev. 1.). We use here focal signals to exemplify the effect of optimal transport. Such signals could correspond to dipolar activations derived from MEG/EEG using dipole fitting methods [START_REF] Scherg | Two bilateral sources of the late AEP as identified by a spatio-temporal dipole model[END_REF] or sparse regression techniques [START_REF] Wipf | Analysis of empirical bayesian methods for neuroelectromagnetic source localization[END_REF][START_REF] Gramfort | Timefrequency mixed-norm estimates: Sparse M/EEG imaging with non-stationary source activations[END_REF].

Figure 1 presents the locations of the two regions (labels), the averages with and without Gaussian smoothing and the Kantorovich average. Gaussian smoothing leads to a highly blurred average which exceeds the extent of the regions of interest, while it also strongly reduces the amplitudes of the signals, potentially washing out the statistical effects. The peak amplitudes obtained with optimal transport are also higher and closer to the individual peak amplitudes. One can clearly observe the limitations of Gaussian smoothing, which furthermore requires to set the bandwidth of the kernel. The Kantorovich average nicely highlights two foci of signals at the group level.

Results on fMRI data. We used here fMRI data analyzed on a voxel grid. It corresponds to 20 subjects from the database described in [START_REF] Pinel | Fast reproducible identification and large-scale databasing of individual functional cognitive networks[END_REF]. We average here the standardized effect of interest induced by left hand button press. In Fig. 2-a we show the Euclidian average without smoothing. In Fig. 2-b we report results obtained by classical averaging following Gaussian smoothing with FWHM of 8 mm. Fig. 2-c shows the Kantorovich mean with constrained mass. One can observe that this barycenter highlights a clear active region without requiring any kernel smoothing. It also leads to a amplitude in the average standardized effect around 1.7 which is much higher than the 0.23 obtained when smoothing.

Results on MEG data. We now evaluate the benefit of the proposed approach on experimental data. These data were acquired with a Neuromag VectorView system (Elekta Oy, Helsinki, Finland) with 306 sensors arranged in 102 triplets, each comprising two orthogonal planar gradiometers and one magnetometer. Subjects are presented with images containing faces of familiar (famous) or unfamiliar persons and so called "scrambled" faces. See [START_REF] Henson | A parametric empirical bayesian framework for the EEG/MEG inverse problem: generative models for multisubject and multimodal integration[END_REF] for more details. Dataset contains 16 subjects. For each one, event related fields (ERF) were obtained by averaging about 200 repetitions of recordings following stimuli presentations. Data were band pass filtered between 1 and 40 Hz. Following standard MEG source localization pipelines [START_REF] Gramfort | MNE software for processing MEG and EEG data[END_REF], a noise covariance was estimated from prestimulus time intervals and used for source reconstruction with the cortically constrained dSPM method [START_REF] Dale | Dynamic statistical parametric neurotechnique mapping: combining fMRI and MEG for high-resolution imaging of cortical activity[END_REF]. The values obtained with dSPM can be considered as F statistics, where high values are located in active regions.

In Figure 3, we present results at a single time point, 190 ms after stimulus onset, which corresponds to the time instant where the dSPM amplitudes are maximum. Data correspond the visual presentation of famous faces. In green, is the border of the primary visual cortex (V1) provided by the FreeSurfer functional atlas. One can observe that the Kantorovich barycenter yields a more focal average nicely positioned in the middle of the calcarine fissure where V1 is located. Such a strong activation in V1 is expected in such an experiment consisting of visual stimuli. To investigate more subtle cognitive effects, such as the response of the fusiform face area (FFA) reported about 170 ms after stimulation in the literature [START_REF] Kanwisher | The fusiform face area: A module in human extrastriate cortex specialized for face perception[END_REF][START_REF] Henson | A parametric empirical bayesian framework for the EEG/MEG inverse problem: generative models for multisubject and multimodal integration[END_REF], we report results obtained on contrasts of ERFs measured after famous faces presentations vs. scrambled faces. As illustrated in Figure 4, Kantorovich mean nicely delineates a focal source of activity in the ventral part of the cortex known as the fusiform gyrus.

These results show that Kantorovich means provides focal activation at the population level despite the challenging problem of inter-subject anatomofunctional variability. They avoid the smearing of the signal or statistical effects of interests which naturally occur when data are spatially smoothed before standard averaging. Note again that here no smoothing parameter with FWHM in millimeters is manually specified. Their solution only depends on the metric naturally derived from the geometry of the cortical surface. With a cortical triangulation containing 10,024 vertices and 16 subjects the computation on a Tesla K40 GPU of one barycenter takes less than 1 min.

Conclusion

The contributions of this work are two-fold. First, by considering non-normalized measures particularly relevant for medical imaging data we extend the current state of the art in barycenter estimation using transportation metrics. Following recent contributions on discrete optimal transport we propose a smoothed version of the transport problem that leads us to an efficient optimization algo-rithm. While many contributions on optimal transport work only in one or two dimensions on a regular grid, our approach can cope with the complex geometry of the brain (irregular grids and surfaces). Only the definition of a ground metric is here required. The algorithm proposed involves simple operations that are particularly adapted to modern GPU hardware and allows us to compute barycenters on full brain data in a few minutes.

Second, with simulations defined on the cortex triangulation, a publicly available fMRI dataset with 20 subjects and an MEG dataset processed with a standard analysis pipeline with 16 subjects we demonstrated the ability of the method to clearly highlight activation foci while avoiding the need to smooth the data. The fMRI data showed a clear activation in the right motor cortex and on the MEG data we showed that the proposed approach better identified activation foci in the primary visual cortex and the fusiform gyrus. Both findings, that are consistent with previous neuroscience literature, show that method proposed yields more accurate results than the current pipelines which furthermore requires to set a kernel bandwidth parameter. The removal of any free parameter in the pipeline is a way towards more reproducible neuroimaging results.

Due to the non-linearity of the approach the estimation of statistical threshold shall be performed with non-parametric permutation tests. When thresholding barycenters as presented in Section 4 it is expected that one will obtain clear clusters.

Notations.

  Let d be the cardinal of Ω. Relabeling arbitrarily all the elements of Ω as {1, • • • , d} we represent the set of non-negative measures on Ω by the non-negative orthant R d + . Let 1 d be the d-dimensional vector of ones. For any vector u ∈ R d we write |u| 1 for the l 1 norm of u, d i=1 |u i |. When u ∈ R d + we also call |u| 1 the (total) mass of vector u. Let M = [m ij ] ∈ R d×d + be the matrix that describes the metric between all locations in Ω, namely m ij = D(i, j), i, j ≤ d. Wasserstein Distance for Normalized Histograms. Consider two vectors a, b ∈ R d + such that |a| 1 = |b| 1 . Both can be interpreted as histograms on Ω of the same mass. A non-trivial example of such normalized data in medical imaging is the discretized ODF used for diffusion imaging data [10]. For p ≥ 1, the p-Wasserstein distance W p (a, b) between a and b is the p th root of the optimum of a linear program, known as a transportation problem [4, §7.2]. A transport problem is a network flow problem on a bipartite graph with cost M p (the pairwise distance matrix M raised element-wise to the power p), and feasible set of flows U (a, b) (known as the transportation polytope of a and b), where U (a, b) is the set of d × d nonnegative matrices such that their row and column marginals are equal to a and b respectively:

1 Ni

 1 entropic regularizer λ > 0, step size c. 2: Compute mean mass ρ = |b j |1. 3: Form virtual cost vector ∆ = quantile(M, q%)1 d . 4: Form augmented d + 1 × d + 1 ground metric M as in Equation (3) 5: Set a = 1 d+1 /(d + 1). 6: while a changes do 7:

Fig. 1 .

 1 Fig. 1. Simulation results with focal random signals generated in areas/labels BA45 (yellow) and MT (red) in a group of 100 subjects. Data are defined on a surface with 10,024 vertices (FreeSurfer fsaverage 5). One shows the standard averaging referred to as Mean, the averaging after Gaussian smoothing is referred to as Mean (S) (mean after Gaussian smoothing with FWHM=8 mm), and the Kantorovich mean (p=1). The result Mean shows the focal signals with random positions in the labels delineated in green. The Kantorovich mean highlights clear foci of activations in the ROIs without smearing the activation as with Gaussian smoothing which furthermore significantly dampens the amplitudes.

Fig. 2 .

 2 Fig. 2. Averaging of the standardized effect of interest on fMRI data. From left to right, a) the Euclidian mean without smoothing, b) the Euclidian mean with smoothing (FWHM=8 mm), c) the Kantorovich mean with Euclidian ground metric and p=1. The later result highlights a clear foci of activations in the ROI without smearing the activation nor damping the amplitudes as much as the kernel smoothing.

Fig. 3 .

 3 Fig.3. Average of dSPM estimates derived from MEG ERF data on a group of 16 subjects stimulated with pictures of famous faces. From left to right: standard mean and Kantorovich mean. The left hemisphere is displayed in medial view. In green is the border of the primary visual cortex (V1) provided by FreeSurfer. One can observe that the Kantorovich mean has its peak amplitude within V1.

Fig. 4 .

 4 Fig. 4. Group averages (16 subjects) of dSPM estimates derived from MEG ERF data obtained by contrasting the famous faces stimulation with the scrambled faces.From left to right: standard mean and the Kantorovich mean. The right hemisphere is displayed in ventral view. Optimal transport results highlight a focal activity in the Fusiform gyrus known to be implicated in face processing[START_REF] Kanwisher | The fusiform face area: A module in human extrastriate cortex specialized for face perception[END_REF].

  |b| 1 and non-empty when |a| 1 = |b| 1 (in which case one can easily check that the matrix ab T /|a| 1 belongs to that set). The p-Wasserstein distance W p (a, b) raised to the power p (written W p p (a, b) below) is equal to the optimum of a parametric Optimal Transport (OT) problem on d 2 variables,

https://surfer.nmr.mgh.harvard.edu/

Acknowledgements A. Gramfort was supported by the ANR grant THALA-MEEG, ANR-14-NEUC-0002-01. M. Cuturi gratefully acknowledges the support of JSPS young researcher A grant 26700002, the gift of a K40 card from NVIDIA and fruitful discussions with K.R. Müller. The work of G. Peyré has been supported by the European Research Council (ERC project SIGMA-Vision).


F Caro

F Coquel

D Jamet

S Kokh

Phase Change Simulation for Isothermal Compressible Two-Phase Flows

We present a numerical scheme based on a two-step convection-relaxation strategy for the simulation of compressible two-phase flows with phase change. The core system used here is a simple isothermal model where stiff source terms account for mass transfer.

I. Introduction

For the past years a great deal of attention has been paid to the simulation of two-phase flows. These flows are involved in a wide range of industrial applications from boiling crisis prediction to cavitating flows. We focus here on the simulation of dynamical phase change between two compressible fluids. Previous works such as Le Metayer 10 managed to achieve such task for the case of very fast mass transfer phenomena. We rather deal here with slower processes as Jamet 8 where both fluids are close to the thermodynamical equilibrium. Moreover, we restrict ourselves to the case of isothermal flows and we also assume that both fluids are in thermal equilibrium.

The phase change model adopted here has been introduced by Caro. 4 This model encompasses two relaxation mechanisms through stiff source terms: a mechanical relaxation effect and a thermodynamical effect that accounts for mass transfer.

We will first briefly recall the basic model and its properties. We shall then describe two equilibrium submodels issued from the relaxations processes when instantaneous equilibrium is assumed. Two numerical schemes will then be introduced. Both schemes are based on a two-step strategy involving a convective step (by the mean of a Godunov-type method) followed by a relaxation step. Finally, we shall present 1D and 2D numerical results with convergence tests.

II. Governing Equations

We consider a flow involving two phases represented by two compressible fluids. Each fluid α = 1, 2 is equipped with an equation of state (EOS) ρ α → f α (ρ α ) verifying

P α = -ρ 2 α df α /dρ α , dP α /dρ α = ρ α dg α /dρ α ,
where ρ α is the density of fluid α, f α , P α and g α = f α + P α /ρ α are respectively the free energy, the pressure and the free enthalpy of fluid α. We note z α , the volume fraction of phase α. We suppose that z 1 + z 2 = 1 and note z = z 1 . Setting m α = z α ρ α , the model we study here reads

         ∂ t m 1 + div(m 2 u) = λ(g 2 -g 1 ), ∂ t m 2 + div(m 2 u)
= -λ(g 2g 1 ), ∂ t (ρu) + div(ρu ⊗ u + P Id) = 0,

∂ t z + u • ∇z = κ(P 1 -P 2 ), (1)
where λ and κ are non negative parameters, ρ = m 1 + m 2 is the global density and P is defined as P = α z α P α . This model shall be referred to as the relaxed system.

A. Hyperbolicity of the Relaxed System

Considering smooth solutions V = (ρ 1 z 1 , ρ 2 z 2 , ρu, z) T , the system (1) is equivalent in 1D to the quasilinear formulation

∂ t V + A(V)∂ x V = R(V), R(V) =      λ(g 2 -g 1 ) -λ(g 2 -g 1 ) 0 κ(P 1 -P 2 )      , A(V) =      uy 2 -uy 1 y 1 0 -uy 2 uy 1 y 2 0 c 2 1 -u 2 c 2 2 -u 2 2u M 0 0 0 u      ,
where y α = ρ α z α /ρ is the mass fraction of the α phase, c 2 α = dP α /dρ α is the squared sound velocity of phase α and M = ∂P /∂z. Denoting by c the mixture sound velocity defined by c 2 = α y α c 2 α , the matrix A(V) possesses three distinct eigenvalues:

λ 1 = u -c, λ 2 = λ 3 = u, λ 4 = u + c. The corresponding left eigenvectors l k and right eigenvectors r k read l 1 = 1 2c 2      c 2 1 + uc c 2 2 + uc -c M      , l 2 = 1 c 2      -y 1 c 2 1 + c 2 -y 1 c 2 2 0 -y 1 M      , l 3 = 1 c 2      -y 2 c 2 1 c 2 -y 2 c 2 2 0 -y 2 M      , l 4 = 1 2c 2      c 2 1 -uc c 2 2 -uc c M      , r 1 =      y 1 y 2 u -c 0      , r 2 =      1 0 u -c 2 1 /M      , r 3 =      0 1 u c 2 2 /M      , r 4 =      y 1 y 2 u + c 0      .
This eigenstructure shows that the system (1) is hyperbolic.

The system (1) is endowed with a (mathematical) entropy-entropy flux pair [ρf + ρ|u| 2 /2, (ρf + ρ|u| 2 /2 + P )u]. Although (ρ 1 z 1 , ρ 2 z 2 , ρu, z) → (ρf + ρ|u| 2 /2) is not convex, the following entropy equation is available

∂ t ρf + ρ |u| 2 2 + div ρf + ρ |u| 2 2 + P u = -λ(g 1 -g 2 ) 2 -κ(P 1 -P 2 ) 2 ,
which allows to derive an entropy inequality associated with the entropy-entropy flux pair The κ = +∞ assumption in the system (1) is formally equivalent to compute z ∈ [0, 1] such that

∂ t ρf + ρ |u| 2 2 + div ρf + ρ |u| 2 2 + P u ≤ 0.
P 1 m 1 z = P 2 m 2 1 -z , (2)
for given fixed values m 1 = ρ 1 z 1 and m 2 = ρ 2 z 2 .

Proposition 1. We suppose that the functions P α : ρ α → P α (ρ α ) are C 1 (R + ), are strictly increasing on R + and tend to +∞ when ρ α → +∞. Then there exists a unique z * (m 1 , m 2 ) ∈ [0, 1] solution of the equation [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF]. Moreover

m 1 → z * (m 1 , m 2 ) is non-decreasing on R + , m 2 → z * (m 1 , m 2 ) is non-increasing on R + , (m 1 , m 2 ) → z * (m 1 , m 2 ) is as regular as ρ α → P α , α = 1, 2.
Using z * , the M-equilibrium system reads

     ∂ t m 1 + div(m 1 u) = λ(g 2 -g 1 ), ∂ t m 2 + div(m 2 u) = -λ(g 2 -g 1 ), ∂ t (ρu) + div(ρu ⊗ u + P Id) = 0, (3)
with P = α z * α P α = P 1 = P 2 . Considering smooth solutions W = (ρ 1 z 1 , ρ 2 z 2 , ρu) T of (3), the system (3) is equivalent in 1D to

∂ t W + B(W)∂ x W = S, B(W) =    u 0 ρ 1 z 1 0 u ρ 2 z 2 (1/ρ)∂ m1 P (1/ρ)∂ m2 P u    , S =    λ(g 2 -g 1 ) -λ(g 2 -g 1 ) 0    .
The matrix B(W) possesses three disctinct eigenvalues:

λ 1 = u -c * , λ 2 = u, λ 3 = u + c * , where c * is defined by 1 ρ(c * ) 2 = α (z * α ) 2 m α c 2 α .
This ensures the strict hyperbolicity of the M-equilibrium system (3).

C. Limit System κ = +∞, λ = +∞ (MT-equilibrium)

The condition κ = +∞, λ = +∞ (MT-equilibrium) is formally equivalent to determine m 1 , m 2 and z, for a given ρ, such that

P 1 m 1 z = P 2 m 2 1 -z , g 1 m 1 z = g 2 m 2 1 -z , ρ = m 1 + m 2 . (4)
System ( 4) is equivalent to find ρ 1 > 0, ρ 2 > 0 and z ∈ (0, 1) such that for a given ρ > 0

P 1 (ρ 1 ) = P 2 (ρ 2 ), g 1 (ρ 1 ) = g 2 (ρ 2 ), ρ = zρ 1 + (1 -z)ρ 2 . ( 5
)
Let us note that the above system does no longer depend on m 1 , m 2 and z, but solely on the parameter ρ. Unfortunately, on the contrary of system (2), nor existence nor unicity for the solution of ( 5) is ensured. However, this issue seems natural. Indeed, it expresses the fact that it is not possible to prescribe a mass transfer mechanism between two fluids with two arbitrary EOS. A certain degree of compatibility between both fluids EOS is required in order to define correctly the thermodynamical equilibrium. This matter can be illustrated in the case of stiffened gas EOS. For such fluids we have where the constants a α and b α are related to the constant flow temperature T as follows

P α = a α ρ α -P ∞ α , g α = b α + a α log(ρ α ), 3
a α (T ) = (γ α -1)c vα T, b α (T ) = q α + γ α c vα -q α + (γ α -1)c vα log (γ α -1)c vα -c vα log T T,
the parameters γ α > 1, P ∞ α ≥ 0, c v α > 0, q α and q α being constant for each fluid α = 1, 2. In this particular case we have the following solvability condition Proposition 2. For given m 1 , m 2 and z, ρ * 1 and ρ * 2 solutions of (5) verify

ρ * 2 = (ρ * 1 ) a1/a2 exp b 1 -b 2 a 2 and a 1 ρ * 1 + (P ∞ 1 -P ∞ 2 ) -a 2 (ρ * 1 ) a1/a2 exp b 1 -b 2 a 2 = 0. ( 6
)
If (a 1 -a 2 )(P ∞ 2 -P ∞ 1
) ≥ 0 then system (6) admits a unique solution, while if (a 1a 2 )(P ∞ 2 -P ∞ 1 ) < 0 it has two solutions or no solution.

Remark 1. For the specific case of perfect gas EOS, system (5) admits a unique explicit solution which reads

ρ * 1 = exp - b 1 -b 2 a 1 -a 2 a 1 a 2 a 2 a 1 -a 2 , ρ * 2 = exp - b 1 -b 2 a 1 -a 2 a 1 a 2 a 1 a 1 -a 2 .
In the sequel we shall consider that the system (5) possesses at most a unique solution ρ * 1 > 0, ρ * 2 > 0 and z * ∈ (0, 1) for a given ρ > 0. Consequently ρ * 1 and ρ * 2 have fixed values which only depend on the fluids EOS definition and, supposing that ρ * 1 < ρ * 2 , z * depends on ρ as follows

z * = 1 if ρ < ρ * 1 , z * = ρ -ρ * 2 ρ * 1 -ρ * 2 if ρ * 1 ≤ ρ ≤ ρ * 2 , z * = 0 if ρ * 2 < ρ. (7)
Finally, the MT-equilibrium system is formally equivalent to

∂ t ρ + div(ρu) = 0, ∂ t ρu + div(ρu ⊗ u + P (ρ)Id) = 0, (8)
where the pressure P is determined by

P (ρ) = P 1 (ρ) if z * = 1, P (ρ) = z * P 1 (ρ * 2 ) + (1 -z * )P 2 (ρ * 2 ) if z * ∈ (0, 1), P (ρ) = P 2 (ρ) if z * = 0.
The pressure law ρ → P (ρ) is continuously differentiable except in points ρ * 1 and ρ * 2 where it is only continuous. System (8) has a structure close to the isothermal Euler equation, however the system is only weakly hyperbolic. Indeed, for ρ / ∈ {ρ * 1 , ρ * 2 } it is possible to compute the Jacobian associated with the system. It possesses two eigenvalues

λ 1 = u + c and λ 2 = u -c, where c 2 = dP/dρ. While c > 0 for ρ / ∈ [ρ * 1 , ρ * 2 ], in the case ρ ∈ (ρ * 1 , ρ *
2 ), we have c = 0. In this case both eigenvalues collapse and the system looses its eigenvector basis.

Remark 2. The MT-equilibrium system is endowed with a volumic specific free energy F which has the following properties: ρ → F (ρ) is only weak-convex, pressure P is related to F by P = ρdF/dρ -F and the graph of

F in the (ρ, F ) is the convex hull of the function ρ → min[F 1 (ρ), F 2 (ρ)],
where F α is the volumic free energy of the fluid α.

Remark 3. The solution of the Riemann problem for the MT-equilibrium system may not be unique in the class of entropic solutions. In this sense, the sytem is ill-posed.

III. Numerical Treatment

From now on, we shall suppose that both fluids are governed by a stiffened gas EOS. Consider i ∈ Z and n ∈ N, let ∆x > 0 and ∆t > 0 be the space and time steps. Recalling that W = (ρ 1 z 1 , ρ 2 z 2 , ρu) T , we note W n i the finite volume approximate value of (1/∆x) (i+1/2)∆x (i-1/2)∆x W(x, t n )dx where t n = n∆t. The numerical strategy we adopt here is a two-step scheme which consists in a first convective step thanks to a discretization of the relaxed system (1) without source terms

W n+1/2 i -W n i + ∆t ∆x F i+1/2 (W n , z n ) -F i-1/2 (W n , z n ) = 0, z n+1/2 i -z n i + ∆t ∆x H i (W n , z n ) = 0, (9)
followed by a projection (or relaxation) step onto the equilibrium states

(W n+1 , z n+1 ) = Π(W n+1/2 , z n+1/2 ).
The convection step ( 9) is achieved in similar lines as Chanteperdrix 5 and Allaire 2 using here a Roe-type linearization for F i+1/2 and an upwind scheme for H i (W n , z n ). For the relaxation step, we propose two variants depending on the chosen equilibrium states set.

A. Relaxation Towards The M-Equilibrium (M-scheme)

For the M-equilibrium, we suppose λ to be a finite but very large parameter. We formally set κ = +∞ which imposes that the equilibrium value z * is retrieved by solving equation [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF]. With the stiffened gas assumption we have

z * = β 1 + β , β = -a 2 m 2 -a 1 m 1 -P ∞ 2 +P ∞ 1 + a 2 m 2 -a 1 m 1 -P ∞ 2 +P ∞ 1 2 +4a 1 a 2 m 1 m 2 2a 2 m 2 . ( 10
)
Using a splitting operator method, we integrate for t ∈ [0, ∆t] the following ODE system

∂ t m 1 = λ g 2 (ρ 2 ) -g 1 (ρ 1 ) , ∂ t ρ = 0, ∂ t (ρu) = 0, (W, z)(0) = (W, z) n+1/2 . ( 11
)
For stability purpose, instead of a direct numerical integration of system [START_REF] Truskinovsky | Kinks versus shocks[END_REF] we seek for an approximate ODE system that allows explicit integration. In order to do so we scale the parameter λ by setting λ = m 1 m 2 λ, where λ is a supposedly large constant parameter. Then, we also freeze the (g 2g 1 ) term in its value (g 2g 1 ) n+1/2 evaluated at (W, z) n+1/2 . We obtain the simplified ODE system

∂ t (m 1 ) = λm 1 (ρ n+1/2 i -m 1 )(g 2 -g 1 ) n+1/2 i , ∂ t ρ = 0, ∂ t (ρu) = 0. ( 12
)
This system can be solved explicitely

m 1 (t) = ρ n+1/2 i ρ n+1/2 i (m 1 ) n+1/2 i -1 exp[ λρ n+1/2 i (g 2 -g 1 ) n+1/2 i t] + 1 -1 , ρ(t) = ρ n+1/2 i and u(t) = u n+1/2 i .
Let us notice that the above expression ensures m 1 (t) > 0. We can now complete the source term integration step by setting

ρ n+1 i = ρ n+1/2 i , u n+1 i = u n+1/2 i , (m 1 ) n+1 i = m 1 (∆t, ), (m 2 ) n+1 i = ρ n+1 i -(m 1 ) n+1 i .
Finally we update z by remapping z onto the M-equilibrium states by using expression [START_REF] Metayer | Modélisation et résolution de la propagation de fronts perméables: application aux fronts d'évaporation et de détonation[END_REF] with W n+1 .

of 10

American Institute of Aeronautics and Astronautics Paper YEAR-NUMBER B. Relaxation Towards The MT-Equilibrium (MT-scheme)

We suppose the equilibrium values ρ * 1 and ρ * 2 known either explicitly for the case of perfect gases (see remark 1), either by solving [START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF] with an iterative algorithm for the case of stiffened gas.

The intermediary solution (W, z) n+1/2 obtained after the convective step ( 9) is projected onto MTequilibrium states using relation [START_REF] Helluy | Finite volume simulation of cavitating flows[END_REF]. This leads to

if ρ n+1/2 i < ρ * 1 z n+1 i = 1, (m 1 ) n+1 i = ρ n+1/2 i , (m 2 ) n+1 i = 0, if ρ * 1 < ρ n+1/2 i < ρ * 2 z n+1 i = ρ n+1/2 i -ρ * 2 ρ * 1 -ρ * 2 , ( m
1 ) n+1 i = ρ * 1 z n+1 i , (m 2 ) n+1 i = ρ * 2 (1 -z n+1 i ), if ρ * 2 < ρ n+1/2 i z n+1 i = 0, (m 1 ) n+1 i = 0, (m 2 ) n+1 i = ρ n+1/2 i , and
ρ n+1 i = ρ n+1/2 i , u n+1 i = u n+1/2 i
.

IV. Numerical Results

A. A One Dimensional Isothermal Phase Change Test

We consider a 1 m long domain containing liquid and vapor initially separated by an interface located at x = 0.5 m. We suppose that both phases are perfect gas (characteristics parameters are given in table 1). The temperature is fixed at T = 100 K. The fluids are supposed to be initially at rest, the densities are chosen so that pressure is constant in the whole domain: ρ = ρ * 1 6.05921 × 10 -4 kg.m -3 at the left of the interface and ρ = ρ * 2 33.21726 × 10 -4 kg.m -3 at the right. We suppose that the left boundary is a piston which moves towards left at constant speed u p = -100 m.s -1 . The right boundary is a reflecting wall. The piston motion will generate an accoustic wave travelling from left to right that will reach the interface. This will perturb the thermodynamical equilbrium of the interface and trigger the system relaxation. A phase change process occurs. left right

C v (J.kg -1 .K -1 ) 1 816 1 040 γ = C p /C v 2.

1.43

Table 1. Fluids parameters for the one dimensional piston test.

We use here the MT-scheme with different space steps. Figure 1 displays the profiles of z, pressure P and velocity u at instants t = 6.66 ms and t = 20 ms. The solid lines represent the approximate solution with 5000 cells, the × symbol lines represent the solution with 1000 cells and the symbol lines represent the solution with 100 cells. This allows to check the good convergence behaviour of the solution when the space step goes to zero. Let us note that some small pressure oscillations appear for the rough mesh, nevertheless their amplitude decrease while refining the space step.

We now turn to the convergence of the solution obtained with the M-scheme for large values of λ to MTscheme appromixate solution. Figure 2 shows the profiles computed with both M-scheme and MT-schemes for the variable z and the pressure in the case of the piston test. We use a λ values ranging from 10 to 10 6 . While for low values of λ the solutions are quite different, we see that as expected for λ → ∞ they tend to the same profile. Le us consider the MT-scheme computed solution V ref h as reference approximate solution. Let V λ h be the solution obtained with the M-scheme. We use here a fixed 500 cells mesh. Figure 3 displays the graph of the function ln( λ)

→ ln V ref h -V λ h L 1 .
We numerically check that for λ → +∞ we have

V λ h -V ref h L 1 = O 1/ λ .
¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¡¢

B. Compression of a Gas Bubble

We consider 1 m long square domain discretized over a 100 × 100 cells mesh. A gas bubble is surrounded by liquid in the center of the domain. The radius of the gas bubble is initially r = 25 cm. The EOS used are stiffened gas type whose coefficients are given in table 2. The temperature is set to T = 600 K. The densities of vapor and liquid are obtained by solving the equation ( 5 c v (J.kg -1 .K -1 ) γ = c p /c v P ∞ (Pa) q (kJ.kg -1 ) q (kJ.kg -1 .K We suppose again the left boundary to be a piston which moves towards left at constant speed u p = -100 m.s -1 . The other boundaries are reflecting wall. The figure 4 shows the z profile obtained with the MT-scheme (on the bottom) and the z profile obtained with the homogeneous system (on the top) at various instants. We notice that with the MT-scheme phase change occurs and the vapor bubble liquefies while for the homogenous off-equilibrium system (1) the bubble is just compressed.

V. Conclusion

We have presented in this paper two variants of a convection-relaxation numerical method for computing phase change phenomena within the framework of a simple two-phase compressible isothermal flows model. The mass transfer effects are treated either by a source term integration step, either by an equilibrium states projection step. Numerical tests have been presented which show the ability of the method to capture solutions involving dynamical phase change processes remaining close to the thermodynamical equilibrium. As the equilibrium system is not strictly hyperbolic, stability and convergence behavior of the approximate solutions are a critical issue. Nevertheless, both space convergence and source term integration in the instantaneous relaxation limit and have been successfully challenged in the numerical tests.

2 of 10 American

 210 Institute of Aeronautics and Astronautics Paper YEAR-NUMBER B. Limit System κ = +∞, λ < +∞ (M-equilibrium)

of 10 American

 10 Institute of Aeronautics and Astronautics Paper YEAR-NUMBER

Figure 1 .

 1 Figure 1. Profiles of z, the pressure of the velocity obtained with the MT-scheme. The solid lines represent the approximate solution with 5000 cells, the × symbol lines represent the solution with 1000 cells and the symbol lines represent the solution with 100 cells.

Figure 2 .

 2 Figure2. Profiles of z and the pressure at t = 20 ms obtained with the MT-scheme (solid lines) and the M-scheme ( symbol) for λ ranging from 10 to 10 6 .

Figure 3 .

 3 Figure 3. Convergence behavior for λ → +∞ of the M-equilibrium scheme. The left and right graphs represent respectively ln( λ) → ln z ref hz λ h L 1 and ln( λ) → ln P ref h -P λ h L 1 at t = 20 ms.

  ), which results the following approximate values ρ = ρ * 1 35.56506 kg.m -3 for the vapor and ρ = ρ * 2 686.24091 kg.m -3 for the liquid.

  ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ ¡¢ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ £ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¡£ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤ ¡¤

	2
	1
	0
	-1
	-2
	-3
	-4
	-5
	2	4	6	8	10	12	14

Table 2 .

 2 Fluids parameters for the two dimensional bubble compression.

	-1 )

of 10 American Institute of Aeronautics and Astronautics Paper YEAR-NUMBER

of

American Institute of Aeronautics and Astronautics Paper YEAR-NUMBER

. z profile obtained with the MT-scheme on the bottom and z profile obtained with the homogeneous system on the top; t is varies from t = 0 s to t = 1.08 ms.