Baptiste Magnier

Arezki Aberkane

Philippe Borianne

Philippe Montesinos

Christophe Jourdan

MULTI-SCALE CREST LINE EXTRACTION BASED ON HALF GAUSSIAN KERNELS

Keywords: Crest lines, multi-scale, steerable filter

Crest line extraction has always been a challenging task in image processing and its applications. It is possible to detect ridges and valleys in images using second order filters. In order to estimate crest lines of variable widths, a multiscale analysis of the image is required. In this paper we propose a new ridge/valley detection method in images based on the difference of rotating Gaussian semi filters adapted in a multi-scale process. Due to the directional filters, we obtain a new ridge/valley anisotropic detector enabling very precise detection of ridge/valley of varied widths. Moreover, as the detector filters compute the two directions of crest lines, even highly bended crest lines are correctly extracted. Numerical comparisons with other oriented Gaussian filters and results on real images validate the interest of this method.

CREST LINE EXTRACTION

Crest lines or curves in an image represent narrow ridges or valleys in the image surface. They may define thin nets inside an image describing roads or rivers in satellite images, blood vessels in medical images or roots in underground images. Therefore, finding these dense and thin structures is an important issue in image processing.

In a grey level image I, ridges are maxima, valleys are the minima [START_REF] Haralick | Ridges and valleys on digital images[END_REF]. Indeed, considering the image surface S, in the cartesian space, the following equation defines all pixel coordinates: S(x, y) = (x, y, I(x, y)) t . Let S T (x, y) be the tangential plane to the surface in all points where: S T (x, y) = (1, 0, I x (x, y)) t , (0, 1, I y (x, y)) t , [START_REF] Haralick | Ridges and valleys on digital images[END_REF] with (I x , I y ) the partial derivatives of I respectively along the x and the y axis. Ridges and valleys are given by the points where the values of S(x, y) are maxima or, respectively, minima in the orthogonal direction to the curve at (x, y).

The literature presents several methods for ridges and valleys detection. In image filtering, convolving I with the second derivative of the Gaussian filter in two dimensions, and then removing the non local maxima in the direction maximizing the variance allows a detection of crest lines [START_REF] Ziou | Line detection using an optimal IIR filter[END_REF]. Also, the Hessian matrix consists of the second partial derivatives of the image after convolving with a Gaussian smoothing kernel. Computation of the orthogonal direction of a curve for a pixel can be tackled by finding the eigenvector that corresponds to the maximum absolute eigenvalue of the Hessian matrix [START_REF] Lindeberg | Edge detection and ridge detection with automatic scale selection[END_REF] [4] [START_REF] Cha | Finding curves in sar ccd images[END_REF] or the Weingarten endomorphism matrix [START_REF] Armande | Thin nets extraction using a multi-scale approach[END_REF] of the considered pixel. However, results obtained by these approaches can present increased false detection rate in noisy images. It is mainly because the high pass filtering used for the second derivative is sensitive to the noise level.

In order to detect oriented contours or ridges, Freeman and Adelson proposed to use oriented Gaussian derivative kernels called steerable filters [START_REF] Freeman | The design and use of steerable filters[END_REF]. Thus, the second order oriented Gaussian becomes a crest line detector as illustrated in Fig. 1(a). As represented in Fig. 1(c), the second partial derivative of anisotropic Gaussian kernels introduced by Perona [START_REF] Perona | Steerable-scalable kernels for edge detection and junction analysis[END_REF] constitutes a good estimator of ridges and valleys (note that the Difference of Oriented Gaussian (DOG) with a standard deviation ratio of 1.6 allows a crest lines estimation because it approximates well a second derivate of Gaussian (Laplacian of Gaussian -LoG- [START_REF] Marr | Theory of edge detection[END_REF] [10]). These oriented kernels characterize the pixels belonging to a contour by the fact that their oriented energy is maximum in the orthogonal direction to their local orientation. Afterward, Jacob and Unser [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF] have developed special families of edges and ridges detectors based on Gaussian derivatives which are derived from the criteria used by Canny to design his famous edge detector [START_REF] Canny | A computational approach to edge detection[END_REF]. In particular, second and forth order Gaussian derivatives are presented for the extraction of crest lines (Fig. 1). These criteria are expressed under an integral form and assembled into a quadratic criterion, whose minimization with a Lagrange multiplier method leads to an eigenvalue calculus. These approaches are able to correctly detect large oriented linear structures like straight lines. For anisotropic filtering methods [START_REF] Perona | Steerable-scalable kernels for edge detection and junction analysis[END_REF], the robustness against noise depends strongly on the two smoothing parameters of the filter (length and width). If these parameters increase, the detection is less sensitive to noise, but thin crest lines are considered as noise, and, as for the conventional crest line detectors, the estimated position of the line moves farther away from the true position. Moreover, as shown in Fig. 3 (b), junctions and corners are not well detected with these different methods [START_REF] Magnier | Ridge and valley junctions extraction[END_REF]. Indeed, the detector value is greater near the corner than on the corner itself [START_REF] Magnier | Perceptual curve extraction[END_REF]. Consequently, the precision of the detected ridge or valley points decreases strongly at corner points and for non straight object contour parts.

Based on difference of two half Gaussian smoothing kernels, a new anisotropic crest line detection method has been developed in [START_REF] Magnier | Perceptual curve extraction[END_REF]. This approach is also able to detect junctions and corners thanks to elongated and oriented filters in two different directions (see Fig. 2). Therefore, smoothing with rotating filters means that the image is smoothed with a bank of rotated half anisotropic Gaussian kernels [START_REF] Montesinos | A new perceptual edge detector in color images[END_REF] [16]:

G (µ,λ) (x, y) = 1 2πλµ • H y e -x 2 2•λ 2 -y 2 2•µ 2 , (2)
where (x, y) are pixel coordinates and (µ, λ) the standarddeviations of the Gaussian filter. Contrary to [START_REF] Perona | Steerable-scalable kernels for edge detection and junction analysis[END_REF], as we need only the causal part of the filter (illustrated on Fig. 1(d) and 2), we simply "cut" the smoothing kernel by the middle, this operation corresponds to the Heaviside function H y along the y axis [START_REF] Montesinos | A new perceptual edge detector in color images[END_REF]. Using a discrete rotation R θ of angle θ of the original image I, by convolution [START_REF] Deriche | Recursively implementing the gaussian and its derivatives[END_REF], we obtain a collection of directional smoothed images

I θ = (R θ • I) * G (µ,λ)
. In a second step, we apply an inverse rotation of the smoothed image and obtain a bank of 360 • /∆θ where ∆θ represents the step between each discretized orientation. This stage is repeated twice with filters of different widths but of the same length before applying the difference: [START_REF] Lindeberg | Scale-space theory in computer vision[END_REF]. Thus, we obtain for each pixel a signal D(x, y, θ) which corresponds to a scan in 360/∆θ directions [START_REF] Magnier | Perceptual curve extraction[END_REF]. Then a ridge/valley operator Σ (µ,λ1) (x, y) is defined as:

D(x, y, θ) = G (µ,λ1) (x, y) - G (µ,λ2) (x, y), with λ 2 /λ 1 = 1.6 [9]
Σ (µ,λ1) (x, y) = D(x, y, θ M1 ) + D(x, y, θ M2 ) + D(x, y, θ m1 ) + D(x, y, θ m2 ) (3)
where θ M1 , θ M2 are the directions of the local maxima of the function D (see Fig. 2) and θ m1 , θ m2 the directions of the local minima. Conditions of detection are as follows: if Σ (µ,λ1) (x, y) > Σ th , the pixel P belongs to a ridge, if Σ (µ,λ1) (x, y) < -Σ th , the pixel P belongs to a valley, where Σ th > 0. Finally, ridges can easily be extracted computing local maxima of Σ (µ,λ1) (x, y) in the η direction, bisector of (θ M1 , θ M2 ) and valleys computing the minima of Σ (µ,λ1) (x, y) in the η direction, bisector of (θ m1 , θ m2 ) (see Fig. 2). Thereafter is developed the method in a multi-scale approach.

1 2 +- - 1 2 _ Y X Y X 1 2
Thin net P M M Fig. 2. Difference of half rotating smoothing kernel description and parameters. It computes the two directions of a curve even if they are discontinuous or highly bended.

MULTI-SCALE THIN NET EXTRACTION

Detection of Crest Lines of Different Widths

Crest lines represent narrow objects in an image, however, these objects happen to have different widths. Thus, secondorder derivatives techniques can be commonly implemented in a multi-scale framework to detect ridges or valleys of different sizes. The choice and number of scales is especially important since too small scales emphasize small details, while producing spurious responses to noise or texture, tangent objects and other features. Otherwise, too large scales deform original objects at detection rounding corners and by loosing details. Scale selection is most often performed by choosing the maximum principal curvature of the image surface among the Gaussian kernels of different scale parameters. According to a function depending on the scale, it is possible to weight the second-order Gaussian filter responses such that this weighted response is maximized when the true scale of the structure is present. In the literature, various measures of crest lines strength have been proposed concerning isotropic Gaussian filters [3] [6]. Oriented Gaussian filters [START_REF] Freeman | The design and use of steerable filters[END_REF] [8] [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF] can be performed in multi-scale framework in order to detect objects of different sizes. Thus, at each pixel of coordinates (x, y), the selected scale W(x, y) is calculated as the maximum response of the oriented operator F θi,σj (x, y), an oriented filter response, at the orientation θ W (x, y) for the discrete sampling of the angle and the scale domain [START_REF] Sánchez | Segmentation of tracheal rings in videobronchoscopy combining geometry and appearance[END_REF]:

     W(x, y) = max j max i F θi,σj (x, y) , θ W (x, y) = arg max θi(x,y) F θi,σj (x, y) , (4)
where θ i (x, y) represents the i th discrete orientation of the filter at a scale σ j maximizing F θi,σj (x, y). As seen in the previous section, σ j = λ 1 for anisotropic DOG filter. Then, crest lines can be extracted computing local maxima of W(x, y) in the direction perpendicular to θ W (x, y). These different multi-scale techniques, isotropic [START_REF] Lindeberg | Edge detection and ridge detection with automatic scale selection[END_REF] [6] [START_REF] Freeman | The design and use of steerable filters[END_REF] [8] [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF] or anisotropic [START_REF] Sánchez | Segmentation of tracheal rings in videobronchoscopy combining geometry and appearance[END_REF] suffer from the problem pointed out in the previous section i.e. they all fail at level of corners.

Based on the analysis of the image singularities, a steerable, multi-scale singularity index for detecting crest lines has been developed [START_REF] Muralidhar | A steerable, multi-scale singularity index[END_REF]. As we will show in the next section, this approach, although robust against light noise, yields many false detections in textured or heavy noisy images.

Multi-Scale Extraction of Crest Lines using Difference of Oriented Half Gaussian Kernels

In order to detect ridges and valleys of different sizes with their associated corners, being robust to noise and preventing false detection in presence of texture, we propose to adapt the last method based on half Gaussian kernels presented in Section 1 [START_REF] Magnier | Perceptual curve extraction[END_REF] in a multi-scale approach. On one hand, the ratio µ/λ 1 must not be so large to enable the detection of bended and faint crest lines (illustrated in Fig. 3 (c -d)). On the other hand, a too short kernel is not robust against noise and does not extract correctly discontinuous lines in the image [START_REF] Magnier | Perceptual curve extraction[END_REF]. For that purpose, fixing µ/λ 1 = 5, preserving the ratio λ 2 /λ 1 = 1.6 (diagrammed in Fig. 4(a)), combining eq. 3 and eq. 4, eventually noting Σ (µ,λ1) (x, y) = F θi,λ1 (x, y), we achieve a multi-scale analysis of the image, computing the maximum response of the anisotropic half Gaussian oriented filter with different scale parameters. Note that in eq. 4, η(x, y) replaces θ i . Thus, from eq. 4 and η(x, y), crest lines can easily be extracted computing local maxima of W(x, y) in the direction η(x, y) followed by a hysteresis operator.

In a pixel precision process, as shown in [START_REF] Armande | Thin nets extraction using a multi-scale approach[END_REF], the optimal scale of the crest line's width called w (i.e. the number of pixels) is λ 1 = w/(2 √ 3). The smallest crest line width corresponds to one pixel, but Shannon's theorem states that the sampling frequency must be greater than twice the highest frequency contained in the signal, so the smallest starting scale corresponds to λ 1 = 0.58. Then, we can fix the scale variation by ∆λ 1 = 1/(2 √ 3) ≈ 0.29.

EVALUATION, RESULTS AND CONCLUSION

Evaluation and Comparison

In order to carry out some quantitative results, we have conducted a number of tests using multi-scale filters with an image of wheat root system including valleys of different widths and sizes. The binary image in Fig. 4 suppressed. The (remaining) connected component is skeletonised, i.e. reduced to single-pixel-wide shapes by a wellknown thinning algorithm [START_REF] Zhang | A fast parallel algorithm for thinning digital patterns[END_REF]. The skeleton is the almost-8-connected-pixel line of local maxima of the distance-tocomponent-border function. The dissimilarity between two pixel lines -skeleton, ridge lines, medial axis, etc -is given by the sum of their differences, i.e. the number of pixels of the first line which are not in the second. The pixel-to-pixel comparison is done in a 3×3 neighborhood centered successively on the position of the first line pixels. For each point (i, j) of the first line must find a point of the second in the 3×3 neighborhood centered on the (i, j) position. The dissimilarity is increased by one each time no point is found. The introduction of neighborhood allows comparing lines with different pixel connectivity and small local variations. The dissimilarity factor is obtained dividing the counted differences by the number of the first line pixels.

In the case of the comparison between a F line to a given R reference line composed from N pixels, the (factor of) dissimilarity allows evaluating the effectives of the true positive, false positive and false negative classes used in the performance [START_REF] Pavez | Objective comparison of contour detection in noisy images[END_REF] or sensitivity [START_REF] Fawcett | An introduction to roc analysis[END_REF] functions which often characterize the produced results. The true positive are defined by the similarity between the two lines i.e. the complement to N of the dissimilarity, the false positive (resp. false negative) by the dissimilarity between R and F (resp. F and R).

With these functions, we analyzed the effect of the several multi-scale oriented crest line detectors on the reference image corrupted by a Gaussian noise, an uniform white noise or adding texture. The uniform noise and the texture are added on the original image (shown in Fig. 4) using the formula: several approaches [START_REF] Freeman | The design and use of steerable filters[END_REF] [8] [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF] adapted in multi-scale (using 19 scales, see eq. 4 and section 2.2), our method and the multi-scale singularity index [START_REF] Muralidhar | A steerable, multi-scale singularity index[END_REF] in function of the texture (i.e. the value of L) or the value of the Gaussian noise standard deviation. Evaluation of the different methods on an image corrupted by an uniform noise is similar with the addition of texture (multi-scale singularity index approach brings a little better result with uniform noise). Steerable Gaussian filters of order 2 or 4 achieved the best results when I m is not much corrupted, however, false negative and false positive points increase if the pollution becomes too heavy. Both of these methods produce interesting results in images containing small noise or texture level because they allow a detection of short crest lines, contrary to anisotropic filters which must have a standard deviation of filter's width five times smaller than its length (i.e. µ/λ 1 = 5). Moreover, steerable anisotropic Gaussian filter shows a lot of false positive points, unlike half Gaussian filters which achieve always the less false negative points detection. Furthermore, the proposed method performs better in the presence of a heavy corruption of I m in terms of false positive point extraction.

I m = (1 -L) • I 0 + L • I N ,

Result on Real Images and Conclusion

In order to validate our approach, we have applied our multiscale filter on real images. Fig. 6(a) represents an image of a river delta captured by a NASA's standard spatial observer, here, we aim to extract all the arms of the river (the wide as the thin). We use 30 scales for each method adapted in multiscale [START_REF] Freeman | The design and use of steerable filters[END_REF] [8] [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF] and ours. Steerable filter methods with full Gaussian kernels are not able to extract all the streams, espe- This method obtains although promising results on large image databases of underground roots and retina medical images. Results on real images validate our method, which is a precise and a robust multi-scale detection approach of ridges and valleys. Next on our agenda is to adapt our method to 3D segmentation of blood vessel [START_REF] Armande | A 3d thin nets extraction method for medical imaging[END_REF] and medical images [START_REF] Prinet | Thin network extraction in 3d images: application to medical angiograms[END_REF].

Fig. 1 .

 1 Fig. 1. Different 2D crest line detectors. For (c) and (d) are used the following ratios: λ 2 /λ 1 = 1.6 and µ/λ 1 = 5.

Fig. 3 .

 3 Fig. 3. Crest line extraction using: (b) full anisotropic Gaussian kernels and (c -d) half Gaussian kernels.

Fig. 4 .

 4 Fig. 4. Addition of texture on a binarized image 600×903.

Fig. 5 .

 5 Fig. 5. Percentage error evaluation of image in Fig. 4(a).

Fig. 6 .

 6 Fig. 6. Multi-scale extraction of river delta.

ACKNOWLEDGEMENTS

The authors would like to thank the French Agence Nationale de la Recherche (project Ecosfix, ANR-2010-STRA-003) and the International Centre for Agricultural Research for Development (CIRAD) for providing the enabling environment and financial resources that allowed us to achieve this work.


Richard Saurel
email: richard.saurel@univ-amu.fr

Carlos Pantano
email: cpantano@illinois.edu

Keywords: diffuse-interface methods, compressible two-phase flow, nonconservative hyperbolic systems, phase transition, interface sharpening

Simulation of compressible flows became a routine activity with the appearance of shock-/contact-capturing methods. These methods can determine all waves, particularly discontinuous ones. However, additional difficulties may appear in two-phase and multimaterial flows due to the abrupt variation of thermodynamic properties across the interfacial region, with discontinuous thermodynamical representations at the interfaces. To overcome this difficulty, researchers have developed augmented systems of governing equations to extend the capturing strategy. These extended systems, reviewed here, are termed diffuse-interface models, because they are designed to compute flow variables correctly in numerically diffused zones surrounding interfaces. In particular, they facilitate coupling the dynamics on both sides of the (diffuse) interfaces and tend to the proper pure fluid-governing equations far from the interfaces. This strategy has become efficient for contact interfaces separating fluids that are governed by different equations of state, in the presence or absence of capillary effects, and with phase change. More sophisticated materials than fluids (e.g., elastic-plastic materials) have been considered as well.

Diffuse-Interface Capturing Methods for Compressible

Two-Phase Flows

INTRODUCTION

In this review, we concentrate on recent advances in the methodology for the simulation of compressible flow with material interfaces. The focus is on advances in the formulation and their efficient numerical implementation in modern computing architectures; finally, we highlight selected recent applications in complex flows. Although we try to be exhaustive, we give precedence to writing a coherent, self-contained description of the state of the art in the field and to discussing certain problems that, in our opinion, deserve immediate consideration. Unintended omissions of contributions to this field are unavoidable, and we apologize beforehand.

In addition to the well-known discontinuous waves that arise in single-phase gas dynamics, a new discontinuity appears as soon as the flow involves two distinct fluids (or materials). This discontinuity is usually named a contact wave and corresponds to the interface that separates the two fluids, e.g., air and water. The interface may have additional thermodynamic properties from those of the two separate phases, which may or may not obey the same type of equation of state (EOS), e.g., ideal gas and stiffened gas (SG). Interface motion has been addressed by several classes of computational methods that can be classified as follows. In the first class, the interface is preserved as a resolved surface, and the grid deforms following the kinematics of the interface. Prominent examples are the Lagrangian [START_REF] Von Neumann | A method for the numerical calculation of hydrodynamic shocks[END_REF] and arbitrary Lagrangian-Eulerian methods [START_REF] Hirt | Arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF]. These methods can be very accurate if the deformation of the interface is not too large, but they become limited by mesh distortions at large deformations. In the second class, the front tracking method [START_REF] Glimm | Three-dimensional front tracking[END_REF]) attempts to reduce these distortions by considering fixed meshes for the bulk flow and tracking the moving interface by Lagrangian markers. The approach requires the management of at least two flow solvers, and it becomes geometrically challenging as the interface distortion grows. In the third class, the development of the level set method [START_REF] Dervieux | A finite element method for the simulation of a Rayleigh-Taylor instability[END_REF], which was adapted to compressible fluids [START_REF] Aslam | Level set methods applied to modeling detonation shock dynamics[END_REF], and the more modern ghost fluid method [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF], used to compute the pressure in mixture cells, allowed for further simplicity and generality. To avoid the complexity related to mesh management with previous methods, the level set method tracks the interface implicitly through an Eulerian function, and two sets of Euler equations were used to store and evolve the fluid variables when needed, in particular in mixture cells. The ghost fluid method is used to transfer the boundary conditions at interfaces from one set of Euler equations to the other. Although apparently simple, this method still requires careful considerations to ensure robustness and conservation of mass. More details can be found in previous reviews [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF][START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF][START_REF] Sethian | Level set methods for fluid interfaces[END_REF].

In this review, we are concerned with computational methods that artificially let the two fluids develop a small amount of mixing at the otherwise resolved interface: This justifies the term diffuseinterface methods (DIMs). This presents its own problems because the thermodynamic state of the mixture needs to be defined as the interface evolves. Typically, the bulk fluids have density and internal energy that are significantly different from that of the mixture, and in general, it is not straightforward to compute the mixture thermodynamics and particularly the pressure across the diffused interface without additional closures. Two subclasses of DIMs are present in the literature. The first considers physically diffuse interfaces, having a well-defined viscocapillary structure [START_REF] Cahn | Free energy of a nonuniform system. I. Interfacial free energy[END_REF]. In order for these methods to work correctly, the spatial resolution must be sufficiently fine to accurately resolve the interface width, i.e., up to a few nanometers. In addition, the phase transition between a liquid and its vapor is described through a cubic-type EOS. In this approach, thermodynamics and capillarity are coupled. To the authors' knowledge, it has never been shown that this approach is capable of modeling immiscible interfaces. This approach is well-suited to the simulation of phase transition in very small systems. The second subclass of DIMs addresses mixture cells having computational origins instead of physical ones. The volume of fluid method [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] in the frame of incompressible fluids was pioneering work in this direction. Here, an extra evolution equation is added to the flow model, representing the volume fraction of a given phase. At this level, the model adopts a two-phase description of the flow, with subvolumes occupied by the phases and several mass balance equations. Saurel & Abgrall (1999a) and [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: reduced equations[END_REF] extended this approach to compressible fluids.

The success of compressible two-phase flow simulations is certainly a consequence of the formidable growth in computational power, readily available to scientists and engineers, combined with the availability of advanced and efficient numerical methods for the simulation of flows possessing discontinuities (a hallmark of compressible flow). These advances have made complex two-phase flow simulations a routine practice in recent years. The great advantage of DIMs is the use of a single set of equations to describe the two-phase nature of the flow at all of the computational cells, avoiding the explicit tracking of the interface separating different phases. This progress has been made possible thanks to the tremendous advances in the development of methods for single-phase compressible flow simulation that took place in the last century. It was realized early that shocks and contact discontinuities could be easily captured (instead of resolved) by the use of enlarged artificial diffusivity without causing errors outside the inner region of the discontinuity that were too damaging (see von [START_REF] Von Neumann | A method for the numerical calculation of hydrodynamic shocks[END_REF]. This is certainly true for simple shock waves, and it allows the prediction of shock propagation in complex situations by using just a few mesh points across the wave inner region. Later, [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF] developed an improved method where cell boundary fluxes were determined using the fundamental Riemann problem structure, which accounts for all types of wave interactions. Many researchers have improved and simplified the wave capturing method [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF][START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF][START_REF] Osher | Upwind difference schemes for hyperbolic systems of conservation laws[END_REF][START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Toro | Restoration of the contact surface in the HLL-Riemann solver[END_REF]. Research is still in active progress to increase accuracy [e.g., essentially nonoscillatory schemes (ENO, WENO) [START_REF] Harten | Uniformly high order accurate essentially nonoscillatory schemes[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF], discontinuous Galerkin methods [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws[END_REF], ADER schemes [START_REF] Titarev | ADER: arbitrary high order Godunov approach[END_REF]] and to improve robustness for strong shocks or contact interfaces [START_REF] Zhang | Positivity-preserving high order finite difference WENO schemes for compressible Euler equations[END_REF].

Contrary to shock-capturing schemes, the computation of interfaces separating two immiscible phases with different thermodynamics has no physical viscous regularization. Consequently, the approaches of von Neumann-Richtmyer- Godunov (von Neumann & Richtmyer 1950[START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF] are not as effective as they are for shocks. Traditional shock-capturing methods applied to contacts generally result in excessive and artificial diffusion. In addition, spurious velocity and pressure fluctuations appear due to the loss of mechanical equilibrium and may result in computational failure. Removing the spurious mechanical phenomena requires a well-behaved relaxation of the mixture thermodynamics as one traverses from one phase to the other (Saurel & Abgrall 1999a[START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: reduced equations[END_REF][START_REF] Lund | A hierarchy of relaxation models for two-phase flow[END_REF]. Furthermore, because the pure phases are governed by hyperbolic systems of equations (supersets of Euler equations), the mixture thermodynamics and governing equations must be compatible with this level of description. Some of these problems have already been addressed (see [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF] for formulations separating two ideal gases), and [START_REF] Shyue | An efficient shock-capturing algorithm for compressible multicomponent problems[END_REF] and Saurel & Abgrall (1999b) considered liquid-gas interfaces and added evolution equations for the SG EOS parameters to compute the termodynamics of the mixture. These methods were generalized and rationalized with the help of multiphase flow modeling (Saurel & Abgrall 1999a[START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: reduced equations[END_REF][START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF][START_REF] Massoni | Some models and Eulerian methods for interface problems between compressible fluids with heat transfer[END_REF][START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF], Saurel et al. 2009). In all of these formulations, the aim is to solve interfaces with a unique set of partial differential equations (an extended flow model) and a unique hyperbolic solver. Furthermore, one retains the philosophy of shock-capturing methods and abandons the tracking or reconstruction of the interface because interface deformations in most compressible flows can become arbitrarily complex.

Finally, it is possible to alleviate the smearing or artificial diffusion introduced by capturing methods using purely numerical techniques that attempt to sharpen the interfaces [START_REF] Shukla | An interface capturing method for the simulation of multi-phase compressible flows[END_REF][START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF][START_REF] Shyue | An Eulerian interface sharpening algorithm for compressible two-phase flow: the algebraic THINC approach[END_REF], Chiapolino et al. 2017c), as well as to increase the order of approximation and global accuracy [START_REF] Loubere | A new family of high order unstructured MOOD and ADER finite volume schemes for multidimensional systems of hyperbolic conservation laws[END_REF]. Additional, physical extensions have been addressed as well, e.g., chemical reactions [START_REF] Petitpas | Modelling detonation waves in condensed energetic materials: multiphase CJ conditions and multidimensional computations[END_REF], phase change [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF], surface tension [START_REF] Perigaud | A compressible flow model with capillary effects[END_REF], solid-fluid interfaces [START_REF] Favrie | Solid-fluid diffuse interface model in cases of extreme deformations[END_REF], and plastic transformation [START_REF] Ndanou | Multi-solid and multi-fluid diffuse interface model: applications to dynamic fracture and fragmentation[END_REF].

The review is organized as follows. In Section 2 we present several flow models of diffuse interfaces and explain the link between them. Section 3 presents the main ingredients for the numerical resolution of some of these models. Section 4 showcases some results to illustrate the methods' capabilities.

THEORETICAL FORMULATIONS

The various DIMs for compressible flows mentioned in Section 1 are based on reductions or enhancements of the nonequilibrium [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF] (BN) model. Here, nonequilibrium refers to the mechanical and/or thermodynamic conditions across the interface; whether it is the former or the latter is made clear below. The BN model is based on the inviscid Euler equations for each pure phase, whereas other sets of conservation laws can be appended to include additional physics. Different models derived from the BN model may change the spatial and temporal scales required to resolve the fields, reduce the number of fields that need to be solved, or enhance the physical description to include other effects not present in the original BN model. Figure 1 highlights the two typical scenarios considered in this review. Figure 1a indicates a two-phase flow composed of a phase that is finely dispersed but otherwise embedded in a carried phase. The typical size of the fine phase is too small for any reasonable computational grid to resolve its fluid a b

Disperse flow Resolved-inferface flow

Typical two-phase flow scenarios addressed by the formulations discussed in this review: (a) disperse flow, represented as a homogenized or averaged mixture, and (b) resolved-interface flow. A prototypical computational grid overlays the two-phase mixture.

and thermochemical details, and therefore, one relies on a homogenized or average description of the flow; this is termed average mixture in the figure. In the second scenario, the size of both phase volumes can be resolved by the available computational mesh, and the interface separating the two phases is accurately captured; this is termed resolved interface in the figure. The conventional notation in two-phase flow literature is used below: α k is the volume fraction, ρ k the density, u k the velocity, p k the pressure, E k = e k + 1 2 u k • u k the total energy, and e k the internal energy, for phase k = 1 or 2. [START_REF] Lund | A hierarchy of relaxation models for two-phase flow[END_REF] provided a rational method to deduce diffuse-interface, two-phase flow models from the total nonequilibrium BN model on the basis of asymptotic analysis. Furthermore, it is possible to model additional physics, e.g., solid mechanics with elastic-plastic flows in similar formulations. We refer readers to [START_REF] Favrie | Solid-fluid diffuse interface model in cases of extreme deformations[END_REF] and [START_REF] Ndanou | Multi-solid and multi-fluid diffuse interface model: applications to dynamic fracture and fragmentation[END_REF] for further details and concentrate below exclusively on two-phase fluid flows.

Subsequently, we will rely on the Noble-Abel stiffened gas (NASG) EOS, because it is convenient and sufficiently general to account for many effects exhibited by compressible fluids:

p k = (γ k -1) ρ k (e k -q k ) 1 -ρ k b k -γ k p ∞,k , 1 .
where γ k , p ∞,k , and b k are characteristic constants of material k, and q k represents the reference energy of phase k; details can be found in Le Metayer & Saurel (2016). This is a generalization of the well-known SG EOS to include covolume effects, i.e., short-distance repulsive effects. The original SG EOS accounts only for thermal agitation through the internal energy contribution and attractive effects through the term -γ k p ∞,k . All basic physical molecular effects are thus considered in the NASG formulation, resulting in an improved thermodynamic description of matter compared to SG. Obviously, other convex EOSs are possible alternatives.

Nonequilibrium Model

In this model, the balance equations of phase k are

∂α k ∂t + u I •∇α k = μ( p k -p k * ), 2. ∂α k ρ k ∂t + ∇ •(α k ρ k u k ) = 0, 3. ∂α k ρ k u k ∂t + ∇ •(α k ρ k u k ⊗ u k ) + ∇(α k p k ) = p I ∇α k + λ(u k * -u k ), 4. ∂α k ρ k E k ∂t + ∇ •[α k (ρ k E k + p k )u k ] = p I u I •∇α k -μp I ( p k -p k * ) + λ(u k * -u k )u I . 5.
From this point onwards, we refer to the right-hand side terms as interaction terms between the phases. Index k * denotes the conjugate phase to k; i.e., k = 1 implies k * = 2, and vice versa. Relaxation parameters μ and λ denote the rate at which pressures and velocities tend to equilibrium, respectively. The quantities p I and u I denote the pressure and velocity of the interfaces, respectively. The set of coefficients, interfacial pressures, and velocities comprises modeling elements, some introduced by [START_REF] Baer | A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials[END_REF], to produce a flow model whose dynamics approximates those of immiscible two-phase flows. Equation 2 represents the k-phase volume fraction α k , transported at velocity u I , defined as

u I = u I + sgn(∇α k ) p k * -p k Z 1 + Z 2 , 6 .
where Z k = ρ k c k is the acoustic impedance and c k is the sound speed of phase k, respectively, and the interfacial velocity existing inside a droplet or bubble cloud is defined by

u I = Z 1 u 1 + Z 2 u 2 Z 1 + Z 2 . 7.
During advection, volume variations due to pressure differences between the phases are considered through the relaxation term μ( p k *p k ), with μ controlling the rate at which pressure equilibrium is reached, defined as

μ = A I Z 1 + Z 2 ,
where A I represents the specific interfacial area of the mixture. Equations 3 and 4 express the mass and momentum balance of phase k, respectively. The velocity relaxation term on the right-hand side of the momentum equation is λ(u k *u k ), where λ is the product of the specific interfacial area with the pressure drag coefficient:

λ = Z 1 Z 2 μ.
Viscous drag effects can be considered as well through correlations based on the Reynolds number. The nonconservative term p I ∇α k represents the pressure force acting at the droplet cloud boundaries. The interfacial pressure at this level is approximated as

p I = p I + Z 1 Z 2 Z 1 + Z 2 sgn(∇α 1 ) • (u 2 -u 1 ), 8.
where the interfacial pressure acting inside the droplet or bubble cloud is defined by

p I = Z 1 p 2 + Z 2 p 1 Z 1 + Z 2 . 9.
The expressions above were obtained by [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF]. Finally, Equation 5 expresses the energy balance of phase k, which is nonconservative due to the term p I u I • ∇α k and the relaxation terms on the right-hand side. These terms appear in the correct form to ensure thermodynamic and mechanical consistency. Equations 6-9 are derived from the local Riemann problem solution that obviously provides correct local interfacial variables, but these linearized solutions should be replaced by nonlinear solutions when dealing with strong shocks.

The system represented by Equations 2-5 is hyperbolic with wave speeds u I , u k , and u k ± c k , where c k denotes the sound speed of phase k. There is one additional wave in comparison with the original BN model because of the interfacial velocity given by Equation 6. Furthermore, each phase admits the following entropy equation:

∂α k ρ k s k ∂t + ∇ • (α k ρ k s k u k ) = 1 T k (Z k + Z k * ) × Z k (Z k + Z k * ) -1 ( p k * -p k ) + sgn(∇α k ) • (u k * -u k )Z k * 2 |∇α k | + μZ k * ( p k * -p k ) 2 + λZ k (u k * -u k ) • (u k * -u k ) . 10.
Note that the right-hand side of Equation 10 is nonnegative for all k; this ensures that the mixture entropy, s = k ρ k α k s k , necessarily agrees with the second law of thermodynamics.

Equations 2-5 are a two-phase model for mixture flows evolving in velocity and temperature nonequilibrium (different temperatures at each side of the interface), e.g., clouds of drops, bubbles, or granular beds; additional physical effects are needed in the latter case [START_REF] Bdzil | Two-phase modeling of deflagration-to-detonation transition in granular materials: a critical examination of modeling issues[END_REF]. To close the system of equations, one needs an estimate for the specific interfacial area A I . The most convenient method is to supplement Equations 2-5 with evolution equations for the specific number density of drops (or grains) n k per unit mass of phase k, obeying the equation

∂α k ρ k n k ∂t + ∇ • (α k ρ k n k u k ) = α k ρ k ṅk ,
where the production rate of the number of drops or bubbles ṅk is modeled by known correlations that depend on the Weber and Reynolds numbers (see [START_REF] Furfaro | Modeling droplet phase change in the presence of a multi-component gas mixture[END_REF]. Knowledge of the number density n k and volume fraction α k enables the calculation of the temporally and spatially dependent drops' radii and associated interfacial areas.

Alternative formulations of the volume fraction transport equation have been proposed recently [START_REF] Saurel | Modelling compressible dense and dilute two-phase flows[END_REF] with the aim of improving the acoustic properties of the BN model, in particular, when one of the phases is discontinuous and incapable of effectively propagating sound (e.g., in dispersed liquid suspensions). For liquid drops suspended in a gas or gas bubbles suspended in a liquid, the volume fraction of the dispersed phase α 1 in Equation 2 is replaced by

∂α 1 ∂t + ∇ • (α 1 u 1 ) = μ( p 1 -p 2 ).
11.

Equation 11 leads to a thermodynamically consistent model, provided that the pressure relaxation coefficient is stiff (μ → +∞). For liquid drops suspended in gas, the wave speeds are u 1 , u 2 , u 2 + c 2 , and u 2c 2 , whereas for gas bubbles suspended in liquid, we have u 2 , u 1 , u 1 + c 1 , and u 1c 1 . Moreover, these alternative formulations remain compatible with the reduced models examined later. In particular, the [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: reduced equations[END_REF] model can be derived from the formulations given by [START_REF] Saurel | Modelling compressible dense and dilute two-phase flows[END_REF].

The above formulation with the described closures is applicable to two-phase mixtures with velocity nonequilibrium, as well as resolved interfaces separating two pure fluids. In the latter situation, volume fractions evolve from 0 to 1 in short distances and the specific interfacial area A I has no meaning. Interface conditions are respected due to the various terms appearing on the right-hand side of the system and the associated interfacial variables. The system is able to fulfill the expected interface condition of mechanical equilibrium (continuity of pressure and normal velocities) under one of two possibilities, λ = μ = 0 or λ = μ → +∞.

In the first possibility, the interface conditions are ensured by the nonconservative terms because the symmetric interfacial variables in Equations 6-8 model contact interface conditions as general solutions of local Riemann problems. This possibility has been used by [START_REF] Layes | Quantitative numerical and experimental studies of the shock accelerated heterogeneous bubbles motion[END_REF] to study shock interaction with a gas bubble. It is also able to deal with permeable interfaces (boundaries of bubbles or droplets clouds). It has been used for permeable granular interfaces by [START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF], who extended the system of Equations 2-5 to account for granular effects. This capability is not possible with the original BN model, where the interfacial variables are estimated as p I = p 2 and u I = u 1 , which turn out to be too crude approximations of the local interfacial pressure and velocity. The second possibility relies on stiff mechanical relaxation. In this limit the mixture evolves with a single pressure and single velocity. This approach was proposed by Saurel & Abgrall (1999a) in a splitting formulation where the hyperbolic system of Equations 2-5 was solved during a time step in the absence of source terms, followed by pressure and velocity relaxation steps with sources μ( p kp k * ) and λ(u k *u k ), respectively, with both λ = μ → +∞ [START_REF] Lallemand | Pressure relaxation procedures for multiphase compressible flows[END_REF][START_REF] Lallemand | Pressure relaxation procedures for multiphase compressible flows[END_REF]. Rather than solving stiff systems of ordinary differential equations, pressure and velocity relaxation solvers are built, resulting in nonlinear algebraic systems free of relaxation parameters. These systems sometimes admit exact solutions when the EOSs of the individual phases are simple enough; otherwise they are solved with a Newton method.

The second possibility based on stiff relaxation solvers is very robust, more so than the first possibility based on nonconservative terms, which has less entropy production. However, the first possibility enables consideration of permeable and impermeable interfaces [START_REF] Saurel | A multiphase model with internal degrees of freedom: application to shock-bubble interaction[END_REF][START_REF] Saurel | Symmetric model of compressible granular mixtures with permeable interfaces[END_REF], as well as velocity disequilibrium two-phase mixtures. When considering situations with resolved interfaces only (the main focus of the present review), as schematized in Figure 1, the method based on stiff relaxation can be simplified to a two-phase mixture model where the phases evolve in both pressure and velocity equilibria.

Reduced Model of Kapila et al.

An asymptotic expansion of the system of Equations 2-5 in the limit of stiff mechanical relaxation (λ = μ → +∞) reduces, to leading order, to the following model due to [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: reduced equations[END_REF]:

∂α k ∂t + u • ∇α k = β k ∇ • u, 12. ∂α k ρ k ∂t + ∇ • (α k ρ k u) = 0, 13. ∂ρu ∂t + ∇ • (ρu ⊗ u) + ∇ p = 0, 14. ∂ρ E ∂t + ∇ • [(ρ E + p)u] = 0, 15.
where the same notation as before is employed, the mixture density is defined by ρ = k α k ρ k , the mixture total energy is given by ρ E = k α k ρ k E k , and the function β k is given by

β k = ρ k * c 2 k * -ρ k c 2 k ρ k c 2 k α -1 k + ρ k * c 2 k * α -1 k * . 16.
When two phases are present, the pressure is given by the mixture EOS:

p(ρ, e, α 1 ) = ρe - α 1 (1-ρ 1 b 1 )γ 1 p ∞,1 γ 1 -1 + α 2 (1-ρ 2 b 2 )γ 2 p ∞,2 γ 2 -1 α 1 (1-ρ 1 b 1 ) γ 1 -1 + α 2 (1-ρ 2 b 2 ) γ 2 -1
. 17.

Here, each phase has been assumed to be governed by the NASG Equation 1, and the mixture Equation 17 is obtained as a consequence of the mixture energy definition and pressure equilibrium p 1 = p 2 . The internal energy e k = e k ( p, ρ k ) is given by the (convex) caloric EOSs as well as the sound speed

c k = c k ( p, ρ k ) of the phases k = 1, 2.
Because the pressure is in equilibrium for this model, each phase is allowed to possess different temperatures T k , which can be determined with the help of the thermal EOS e k (ρ k , T k ), which is not necessary to integrate the Kapila model. Furthermore, the right-hand side term of Equation 12 is precisely the limit of the relaxation term μ( p 1p 2 ) in Equation 2, and the system of Equations 12-15 admits the following entropy equations,

∂α k ρ k s k ∂t + ∇ • (α k ρ k s k u) = 0, 18.
which are thermodynamically consistent. Also, the system of Equations 12-15 is hyperbolic with wave speeds u and u ± c w , where c w denotes the mechanical equilibrium sound speed [START_REF] Wood | A Textbook of Sound[END_REF] given by

1 ρc 2 w = α 1 ρ 1 c 2 1 + α 2 ρ 2 c 2 2 . 19.
Equations 12-15 describe mixtures of compressible fluids in mechanical equilibrium but out of thermal equilibrium. Although the equations were derived to model granular mixtures of reactive materials, they are suitable for the computation of interfaces separating pure fluids. The contact interface conditions of this system have equal pressures and equal normal velocities. In addition, because it involves two temperatures, the absence of temperature relaxation in numerically diffuse interfaces produces pressure oscillation-free solutions. This property is absent in the reactive Euler equations, widely used in combustion, which assume temperature equilibrium among the constituents and tend to produce spurious (artificial) pressure oscillations [START_REF] Abgrall | How to prevent pressure oscillations in multicomponent flow calculations: a quasi conservative approach[END_REF].

Pressure Nonequilibrium Model

The Kapila model introduces numerical difficulties due to the nonconservative term on the righthand side of Equation 12. This difficulty is overcome by the homokinetic model in pressure nonequilibrium (PNE) derived by [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF]. It consists of a hyperbolic overdetermined system of seven partial differential equations, given by

∂α k ∂t + u • ∇α k = μ( p k -p k * ), 20
.

∂α k ρ k ∂t + ∇ • (α k ρ k u) = 0, 21. ∂ρu ∂t + ∇ • (ρu ⊗ u) + ∇ p = 0, 22. ∂α k ρ k e k ∂t + ∇ • (α k ρ k e k u) + α k p k ∇ • u = -μp I ( p k -p k * ), 23. ∂ρ E ∂t + ∇ • [(ρ E + p)u] = 0, 24.
where

p = k α k p k . 25.
The interfacial pressure is defined by Equation 9, and the entropy equations are given by

∂α k ρ k s k ∂t + ∇ • (α k ρ k s k u) = μ( p k -p k * ) 2 Z k Z k + Z k * . 26.
In this formulation, the volume fraction corresponds to a simple transport equation with relaxation, for which there is no difficulty in preserving volume fraction positivity. Another advantage of this model is that it exhibits a mixture sound speed c f defined as

ρc 2 f = k ρ k c 2 k ,
which is monotonic with respect to the volume fraction, contrary to the Kapila model. This monotonic behavior simplifies the construction of numerical methods and better respects wave propagation in diffuse interfaces.

The procedure to solve the overdetermined system of equations is as follows. The internal energies are just used to determine the volume fractions at the end of the pressure relaxation process, and the presence of nonconservative terms in these equations is of minor importance. The latter effects are minor at pressure equilibrium. The volume fraction is then used in the mixture Equation 17 that takes the mixture internal energy as argument, and not that of the phases. The mixture internal energy is deduced from the (conservative) total energy equation. Then, from this mixture pressure p, the phase's internal energies are reset with their equations of state e k = e k ( p, ρ k ), removing a possible skid effect related to nonconservative terms. Thanks to these successive corrections, numerical computations even in extreme conditions have shown excellent convergence to exact solutions for interface separating (nearly pure) fluids.

Additional Physics

The previous formulations can be enhanced to include several important physical effects, including surface tension and phase change, as described below.

2.4.1. Surface tension. Following Brackbill et al. (1992), the surface tension can be incorporated as a body force, defined by

F σ = -σ κ ∇Y k [Y k ] ,
where Y k = (α k ρ k )/ρ is the mass fraction of one of the phases (liquid or gas), [Y k ] is the mass fraction jump across the interface, κ is (twice the mean of ) the curvature of the interface, and σ is the surface tension coefficient, which is assumed constant here. Thanks to this formulation of capillary effects, the Kapila model can be extended [START_REF] Perigaud | A compressible flow model with capillary effects[END_REF]) by adding to the momentum and energy equations the capillary force and power, respectively, and by attaching one equation for the mass fraction of one of the phases, giving

∂α k ∂t + u • ∇α k = β k ∇ • u, 27. ∂ρ ∂t + ∇ • (ρu) = 0, 28. ∂ρY j ∂t + ∇ • (ρY j u) = 0, 29. ∂ρu ∂t + ∇ • (ρu ⊗ u) + ∇ p = -σ κ∇Y j , 30
.

∂ρ E ∂t + ∇ • [(ρ E + p)u] = -σ κu • ∇Y j , 31.
where j = 1 or 2 (but not both because k Y k = 1) and κ = ∇ • (∇Y j / ∇Y j ). Equations 30 and 31 can also be rewritten in conservation form:

∂ρu ∂t + ∇ • (ρu ⊗ u) + ∇ p = ∇ • (σ W), 32. ∂ρ Ẽ ∂t + ∇ • [(ρ Ẽ + p)u] = ∇ • (σ W • u), 33.
where

W = ∇Y j I - ∇Y j ⊗ ∇Y j ∇Y j , 34
.

ρ Ẽ = ρ E + σ ∇Y j . 35.
The total energy E is still given by the definition used in the Kapila model.

It is worth remarking that there are some differences between the present formulation and the second gradient formulation of [START_REF] Cahn | Free energy of a nonuniform system. I. Interfacial free energy[END_REF]. First, the current thermodynamics is independent of capillary effects, by construction. Second, the conservative formulation in momentum and energy guarantees correct pressure jumps at interfaces, independently of mesh resolution.

Phase change.

Phase transition can be considered as well by adding Gibbs free energy relaxation effects to the mass and volume fraction equations:

∂α k ∂t + u • ∇α k = β k ∇ • u + ρ ρ I ν(g k * -g k ),
36.

∂ρY j ∂t + ∇ • (ρY j u) = ρν(g j * -g j ),
37.

where g k = h k -T k s k is the Gibbs free energy of phase k, and ν is the evaporation/condensation kinetic rate. To close the system of equations, one uses the appropriate convex EOS (Equation 1) for each phase. [START_REF] Saurel | Modelling phase transition in metastable liquids: application to cavitating and flashing flows[END_REF] specified the interfacial density ρ I that appears in the volume fraction equation to be

ρ I = ρ 1 c 2 1 α 1 + ρ 2 c 2 2 α 2 c 2 1 α 1 + c 2 2 α 2 -1
.

The specific interfacial area must be specified in two-phase mixtures (e.g., bubbly and droplet flows). The parameter ν is determined by a subscale model based on local saturation conditions and heat and mass diffusion; readers are referred to [START_REF] Abramzon | Droplet vaporization model for spray combustion calculations[END_REF], [START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF]Furfaro &Saurel (2016) for details. The local saturation condition p = p sat (T ), a consequence of g 1 (T, p) = g 2 (T, p), is anyway valid in almost all liquid-vapor interface situations. Therefore, when dealing with resolved interfaces, local thermodynamic equilibrium is enforced by taking ν → ∞. This treatment is simple and accurate enough to determine the amount of liquid and vapor at saturation when the energy needed for phase change is already present in the liquid. This is the case of cavitating and flashing flows, as discussed by [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF]. A method to deal with stiff phase transition (ν → ∞) is summarized in Section 3.

Boiling flows.

When dealing with boiling flows, modeling of additional physics is needed because the energy responsible for phase change is not present initially in the liquid but comes from external boundaries. Adding heat diffusion and buoyancy renders the model appropriate for the computation of boiling flows. However, when heat diffusion is considered, temperature becomes continuous at interfaces, unlike interfaces with simple mechanical contact. The previous formulations involving two temperatures may be reduced to a single temperature model, which is easier to solve because the nonconservative volume fraction equation vanishes. The formulation uses Equation 37 and considers gravity and heat conduction effects in the momentum and energy equations (Le [START_REF] Martelot | Towards the direct numerical simulation of nucleate boiling flows[END_REF], resulting in

∂ρu ∂t + ∇ • (ρu ⊗ u) + ∇ p -∇ • (σ W) = ρg, 38. ∂ρ Ẽ ∂t + ∇ • [(ρ Ẽ + p)u] -∇ • (σ W • u + λ c ∇T ) = ρg • u, 39.
where g is the gravity acceleration vector, λ c = α 1 λ 1 + α 2 λ 2 represents the mixture thermal conductivity, and Ẽ is defined in Equation 35. This system is closed by the mixture EOS

p(ρ, e, Y 1 ) = A 1 + A 2 -(P ∞,1 + P ∞,2 ) 2 + 1 4 A 2 -A 1 -(P ∞,2 -P ∞,1 ) 2 + A 1 A 2 40.
where

A k = Y k (γ k -1) C v,k Y 1 C v,1 + Y 2 C v,2 ρ(e -q ) -P ∞,k
and q = Y 1 q 1 + Y 2 q 2 . This EOS is derived from mixture density and internal energy definitions under the constraints of equal pressures and equal temperatures. One can verify that Equations 36 and 37 are hyperbolic and possess the entropy equation,

∂ρs ∂t + ∇ • ρs u -λ c 1 T ∇T = ρν(g 2 -g 1 ) 2 T + λ c (∇T ) 2 T 2 , 41.
which satisfies the second law of thermodynamics.

NUMERICAL METHODS

The mathematical structure of the previous formulations can be described according to the following generic representation:

∂Q ∂t + ∇ • F(Q) + B(Q) • ∇Q = S(Q), 42.
where Q is the vector of evolution variables (conserved and not conserved), F is a flux function, B is the advection matrix of all differential terms that cannot be written in conservation form (i.e., not in F), and S is a source. In some of the previous formulations, B or S may be identical to zero. It is assumed here that all intermediate quantities-u I , u I , p I , p I , λ, μ, etc.-are expressed as functions of Q and eliminated from the list of dependent variables in order for us to retain only the evolution fields.

The most accepted method to solve Equation 42 for compressible two-phase flows is the finite volume method (FVM). The FVM has proven to be a robust and reliable numerical approach for the simulation of compressible flows in the presence of shock wave interactions. The basis of the FVM is discussed in detail in standard references [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF]. The starting point is a decomposition of the computational domain into disjoint cells C i , which may be associated with a structured or unstructured grid, and one then integrates Equation 42 over each cell and the interval of time advancement, t = t n+1t n , where t n denotes the times at which the solution is calculated, resulting in

Qn+1 i = Qn i - 1 V (C i ) t n+1 t n ∂C i F(Q) • n dSdt - 1 V (C i ) t n+1 t n C i B(Q) • ∇Q dxdt + 1 V (C i ) t n+1 t n C i S(Q) dxdt, 43.
where

Qn i = 1 V (C i ) C i Q(x, t n ) dx, 44.
with V (C i ) denoting the volume and ∂C i the surface of cell C i , respectively. The solution strategy usually employed to completely define Equation 43 uses operator-splitting techniques [START_REF] Glowinski | Splitting Methods in Communication and Imaging[END_REF] according to which the flux, nonconservative, and source integrals are determined sequentially, independently of each other (care must be used if time accuracy higher than first order is desired). The flux term, second on the right-hand side of Equation 43, can be obtained by the Godunov (flux difference-splitting) or the Boltzmann (flux vector-splitting) methods. The details of each of these approaches are beyond the scope of this review; they can be found in standard books on the FVM [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF].

The general strategy is typically as follows: (a) Use Qn i to reconstruct Q(x, t n ) in each cell. (b) Solve exact or approximate Riemann problems at the faces of the computational cells (ignore sources). (c) Determine fluxes and nonconservative product contributions, perform the integrals in Equation 43, and advance to Q * i by a Godunov-type method (or higher-order variant). (d ) Integrate the source contribution (when present) from

Q * i to Qn+1 i , typically assuming 1 V (C i ) t n+1 t n C i S(Q) dxdt ≈ t n+1
t n S( Qi (t)) dt. 45.

A difficult challenge in the simulation of compressible two-phase flows in the presence of shocks or discontinuities is the nonconservative term, the third term on the right-hand side of Equation 43. For Godunov-type methods based on Riemann problem solutions, the structure of the solution is determined by both F and B; i.e., operator splitting cannot be used safely. The nonconservative structure precludes the unique determination of wave speeds, and additional kinematic conditions must be attached to the formulation (see Dal [START_REF] Dal Maso | Definition and weak stability of nonconservative products[END_REF]. Only when shocks are not present can one resort to the operator-splitting method and discretize the nonconservative integral directly. In the next section, we describe the general computational aspects that are pertinent to the different models.

Methods for Baer-Nunziato-Type Models

BN-type models are perhaps the most delicate systems to solve among those presented above. They are of particular interest when velocity and temperature are different at each side of the interface. The Riemann problem has up to six or seven wave speeds depending on the variants of the formulation, rendering its solution intricate. In addition, many nonconservative terms are present, and the evolution of the nonconservative volume fraction cell average has no physical sense in the presence of discontinuities. Several approaches have been developed and improved over the years to address these difficulties. Saurel & Abgrall (1999a) constructed a Godunov-type method for these equations, in conjunction with a treatment of nonconservative terms, that preserves mechanical equilibrium at interfaces. They developed stiff pressure and velocity relaxation solvers to ensure interface conditions of equal pressures and normal velocities, allowing the model to deal with diffuse interfaces. [START_REF] Abgrall | Discrete equations for physical and numerical compressible multiphase mixtures[END_REF] introduced two improvements. First, the seven waves present in the flow model are considered through a composite two-phase Riemann solver, consisting of several Riemann solvers for the Euler equations. Second, the nonconservative terms are solved accurately by enforcing that a local diffraction occurs each time a nonlinear wave (shock or expansion) interacts with an interface, rendering both pressure and velocity locally constant precisely at the location where the volume fraction is discontinuous. Thanks to this diffraction process, computed with the help of a Riemann solver for the Euler equations, the various nonconservative products are computed unambiguously. The only remaining issue is the interpretation of the volume fraction cell average at shocks. Saurel et al. (2007a) and [START_REF] Petitpas | A relaxation-projection method for compressible flows. Part II: artificial heat exchanges for multiphase shocks[END_REF] studied this problem further but a general solution is still lacking. [START_REF] Tokareva | HLLC-type Riemann solver for the Baer-Nunziato equations of compressible two-phase flow[END_REF] and [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF] have proposed methods that use the Harten-Lax-van Leer-contact (HLLC)-type Riemann solver including nonconservative products, whereas [START_REF] Dumbser | FORCE schemes on unstructured meshes II: non-conservative hyperbolic systems[END_REF] and [START_REF] Franquet | Runge-Kutta discontinous Galerkin method for the approximation of Baer and Nunziato type multiphase models[END_REF] developed discontinuous Galerkin methods for the BN model. The methods that follow consider flow models with a single interface velocity.

Methods for the Kapila Model and Its Variants

The flow model of Kapila is particularly interesting for simulating interfaces separating fluids.

It also has the ability to create interfaces dynamically. This capability is a consequence of the nonconservative term on the right-hand side of Equation 12. When positive, the velocity divergence term results in subscale bubble growth, giving the appearance of macroscopic interfaces (e.g., gas cavities). However, this model has two important numerical difficulties. First, under some conditions, the nonconservative divergence terms just mentioned may produce negative volume fractions, resulting in immediate simulation failure. Second, the manner by which shock energy is partitioned between the phases may be ambiguous. The PNE formulation of [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] addressed the first issue by replacing the nonconservative term in Equation 12with the relaxation one in Equation 20. This was achieved at the cost of inserting extra energy equations to compute pressure imbalance. The method of [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] fits in the frame of Equation 43, where a specific HLLC-type Riemann solver is built based on the sound speed of the PNE model. The hyperbolic step is followed by pressure relaxation to reach the equilibrium pressure of the mixture. The method is simple and very robust, even in extreme flow conditions.

A second issue for the Kapila model is related to shock relations. These relations are needed when a shock propagates in a two-phase mixture, not when the interface separates two pure (or nearly pure) fluids. In the single-phase limit, the correct shock speeds are computed, because the model tends to the Euler equations. Difficulties arise when the two phases are in nonnegligible proportions (e.g., solid alloys, composite materials, or granular mixtures). Saurel et al. (2007b) obtained these relations in the weak shock limit. The corresponding system of jump relations reads

Y k = Y 0 k , 46. ρ(u -) = ρ 0 (u 0 -) = m, 47. p -p 0 + m 2 (v -v 0 ) = 0, 48. e k -e 0 k + p + p 0 2 (v k -v 0 k ) = 0, 49.
where is the shock speed and v = 1/ρ the specific volume. Variables with superscript 0 denote the state upstream of the shock. A Riemann solver based on Equations 46-49 can be found in [START_REF] Petitpas | A relaxation-projection method for compressible flows. Part II: artificial heat exchanges for multiphase shocks[END_REF], but imposing the shocked state in a given unsteady computation is challenging. Indeed, the cell average of the volume fraction produces skid effects and the (strong) shock becomes a succession of weak shocks. [START_REF] Saurel | Simple and efficient relaxation methods for interfaces separating compressible fluids, cavitating flows and shocks in multiphase mixtures[END_REF] derived a method based on interphase artificial heat exchanges to correctly partition the shock energy among the phases. [START_REF] Petitpas | Modelling detonation waves in condensed energetic materials: multiphase CJ conditions and multidimensional computations[END_REF] derived another simpler, more general method that relies on extra transport equations relative to the Hugoniot pole ( p 0 , v 0 k ), given by

∂ p 0 ∂t + u • ∇ p 0 = 0, ∂v 0 k ∂t + u • ∇v 0 k = 0. 50.
These pole coordinates are used in the energy jump relation e ke 0 k + ( p + p 0 )(v k -v 0 k )/2 = 0 and are directly used in the pressure relaxation solver, based on

e k ( p, v k ) -e 0 k ( p 0 k , v 0 k ) + p(v k -v 0 k ) = 0, k = 1, 2,
with α 1 + α 2 = 1, where p = p (the relaxed pressure) at any point and p = ( p + p 0 )/2 inside the shock layer that needs detection with the PNE flow model. [START_REF] Petitpas | Modelling detonation waves in condensed energetic materials: multiphase CJ conditions and multidimensional computations[END_REF] and [START_REF] Schoch | Multi-phase simulation of ammonium nitrate emulsion detonations[END_REF] give the various details. However, this method has an important limitation because the considered shock must propagate in a uniform mixture. Imposing given shock conditions in this nonconservative flow model with general initial conditions is still an open research area. But it is not a limitation for the computation of flows with interfaces between (nearly) pure fluids.

Methods for Models in Both Temperature and Pressure Equilibria

These models are used in cavitating and boiling flows. The additional temperature equilibrium condition implies that one of the partial differential equations is not needed. The best option is to remove the volume fraction equation that is nonconservative. The flow model now involves four partial differential equations (two mass equations, one mixture momentum equation, and one mixture energy equation) and is fully conservative and hyperbolic. Excluding surface tension terms that need special care [START_REF] Perigaud | A compressible flow model with capillary effects[END_REF], the numerical approximation of this flow model is straightforward. The HLLC Riemann solver is recommended and is used with the speed waves of the Kapila model, which are easier to compute than that of the temperature-pressure equilibrium model. Better stability and robustness are obtained because these wave speeds are greater than that of the temperature-pressure equilibrium model. Details of the overall algorithm are provided by [START_REF] Saurel | A general formulation for cavitating, boiling and evaporating flows[END_REF].

Approach for Phase Transition

For given mass and energy of a two-phase mixture in both temperature and pressure equilibrium, the requirement that the Gibbs free energy is equal for the two phases,

g 1 (T , p) = g 2 (T, p), 51.
provides enough information to determine the thermodynamic state and concentrations. The specific volume and internal energy of the mixture can be expressed in terms of the saturation pressure p = p sat (T ) or temperature T = T sat ( p), given by

v mix = Y 1 v sat 1 ( p) + (1 -Y 1 )v sat 2 ( p),
52.

e mix = Y 1 e sat 1 ( p) + (1 -Y 1 )e sat 2 ( p). 53.
This system involves two unknowns, p and Y 1 , and it is simple to solve unless the thermodynamic state crosses the phase diagram boundaries into one of the single phases. Rather than solving this nonlinear system, a particularly simple and robust alternative method has been proposed by Chiapolino et al. (2017b). It is based on a minmod-type procedure applied to the mass transfer source term of Equation 37 expressed equivalently as

∂ρY 1 ∂t + ∇ • (ρY 1 u) = ρν(Y * 1 -Y 1 ), 54.
where Y * 1 is the mass fraction at thermodynamic equilibrium and ν is a relaxation parameter that is taken to be very large. From the definitions given in Equations 52 and 53, the following estimates of the mass fraction at equilibrium are determined:

Y m 1 = v mix -v sat 2 v sat 1 -v sat 2 , Y e 1 = e mix -e sat 2 e sat 1 -e sat 2 .

55.

The variations between the current local mass fraction are computed as

δY m 1 = Y m 1 -Y 1 and δY e 1 = Y e 1 -Y 1 .
These two variations may have the same sign or opposite sign, and their product, χ = δY m 1 δY e 1 , can be used to estimate the mass fraction at equilibrium by

Y * 1 = Y 1 , i f χ < 0, Y 1 + sgn(δY m
1 ) min( δY m 1 , δY e 1 ), otherwise.

56.

The mass fraction Y 1 is then set to the equilibrium estimate value Y * 1 . This method is particularly simple to manage, even at single phase boundaries. An extension to multicomponent gas mixtures is given by Chiapolino et al. (2017a).

Minimizing Numerical Smearing

DIMs based on the models presented above are simple and efficient in a wide range of complex flow situations, even at high-pressure, -temperature, and -velocity conditions (see Section 4). However, their numerical implementation with shock-/interface-capturing numerical methods tends to produce solutions that exhibit excessive numerical diffusion. This can be seen clearly in low-speed transient conditions. Several approaches have been developed to alleviate this weakness. Here, we highlight some of the solutions that have been successfully used to prevent the artificial smearing of contact interfaces by capturing methods in resolved-interface simulations. For simplicity, we refer to α k as φ below, unless stated otherwise, and consider the behavior of the volume fraction-governing equation surrounding an interface, which is usually detected by requiring that

φ(1 -φ) > ,
where is typically chosen in the range 10 -2 -10 -4 .

Interface compression technique.

The idea behind the compression technique is to modify the volume fraction equation (here for the Euler model) as follows:

∂φ ∂t + u • ∇φ = U o n • ∇ [ h |∇φ| -φ(1 -φ)] , 57.
where U o is a mobility parameter, usually taken to be large with respect to the characteristic speeds in the flow, h is a length scale that sets the desired thickness of the interface, whose chosen value should be larger than the grid spacing, and n = ∇φ/|∇φ| is the local normal vector of the interface. The term on the right-hand side of Equation 57 dominates when U o is large and drives the profiles of φ to the equilibrium profile, defined by the vanishing of the term in square brackets in Equation 57, given by

φ = 1 2 1 + tanh n 2 h ,
where n denotes a normal coordinate to the interface. This idea was first proposed by [START_REF] Olsson | A conservative level set method for two phase flow II[END_REF] in the context of incompressible flow, extended by [START_REF] Shukla | An interface capturing method for the simulation of multi-phase compressible flows[END_REF] to compressible flow, and improved for thermodynamic consistency and extended to the Kapila model by [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF]. Figure 2a,d shows a comparison of a shock in water, generated by a region of 40-MPa overpressure, interacting with an air cavity without and with interface compression, respectively [START_REF] Shukla | An interface capturing method for the simulation of multi-phase compressible flows[END_REF]. Improvements in the solution with interface compression are self-explanatory. In this example, the compression terms were discretized with a second-order centered scheme. [START_REF] Shukla | Nonlinear preconditioning for efficient and accurate interface capturing in simulation of multicomponent compressible flows[END_REF] (c,f ) Unstructured grid simulation results of shock interaction with a gas bubble. In each row, the same resolution was used in corresponding left and right simulations, with the shock traveling from left to right, but unfortunately the conditions and resolution varied from the top row to the bottom row.

Xiao 2014). This technique was first applied to the [START_REF] Allaire | A five-equation model for the simulation of interfaces between compressible fluids[END_REF] model. The key idea of this method is to replace the linear reconstruction of φ(x) from the cell-averaged φi , used in secondorder MUSCL-type (monotonic upstream-centered scheme for conservation laws) methods, with a nonlinear reconstruction of the form

φ(x) = 1 2 1 + tanh β x -x i-1/2 x -xo , 58.
here exemplified for a one-dimensional (1D) arrangement in the cell [x i-1/2 , x i+1/2 ] with φi increasing in the positive x direction, where xo is the unknown location of the interface and β ≈ 2 is a user parameter. Imposing the conservation constraint on Equation 58, given by φi First-order Overbee (solid line) and second-order Superbee (dashed line) limiters and the boundaries of the first-and second-order total variation diminishing (TVD) regions.

= 1 x x i+1/2 x i-1/2 φ(x) dx, 59.
on the cell C i with x = x i+1/2x i-1/2 and inverting for the unknown xo , we have xo = tanh -1 α -δ exp(β(2 φi -1)) , 60.

where α = (tanh β) -1 and δ = (sinh β) -1 . Once xo is known, φ(x) can be evaluated at any x from Equation 58, as required to determine numerical fluxes. This carries over transparently to highorder methods. An example of the method applied to a Mach 1.22 shock wave traveling over an R22 gas bubble is shown in Figure 2b,e.

3.5.3. Total variation diminishing limiter technique. Chiapolino et al. (2017c) described a third method recently developed to sharpen interfaces in unstructured meshes. This approach can be extended to an arbitrary number of fluid components and is based on the MUSCL reconstruction scheme applied to the system of Equations 20-24. First, the limiter function of all fields is set to zero in the interfacial zone, and second, a first-order total variation diminishing (TVD) limiter is used for φ. Figure 3 shows the classical representation of the limiter function θ as a function of the slope ratio r (see [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] for discussion about limiters). This approach is unusual in the context of MUSCL-type schemes, where gradient limiters always belong to the second-order TVD region. Here, the limiter used for volume fraction computation partly resides outside the second-order region and consequently is inappropriate for smooth flows, but behaves very well for discontinuous fields like the volume fraction at interfaces. Figure 2c,f shows results from an example simulation for a Mach 1.5 shock interacting with a Krypton air bubble using an unstructured triangle-based mesh. The interesting point to note is the difference between computed results with the Superbee limiter (Figure 2c) and the new Overbee limiter (Figure 2f ), the latter showing clear improvements. When time evolves, the solution with the Superbee limiter becomes increasingly diffused, whereas that obtained with the Overbee limiter remains sharp.

a b

Low-order method High-order method

Volume fraction iso-surface comparison of collapsing bubble with (a) low-order (i.e., second-order HLLC) and (b) high-order (i.e., fifth-order, weighted essentially nonoscillatory) methods. The same grid was used in both simulations, and the volumes of the two bubbles are nearly identical. Adapted with permission from [START_REF] Tiwari | A diffuse interface model with immiscibility preservation[END_REF]. Abbreviation: HLLC, Harten-Lax-van Leer-contact.

Numerical Anisotropy

The methods described in previous sections can be used in conjunction with structured or unstructured grids and for problems concerned with mixtures or resolved interfaces. In addition, one can use the formulations previously discussed for simulations of coarse (engineering) or fine (research) levels of detail or complexity. If one is interested in highly accurate simulations, especially for resolved interfaces, it is desirable to employ advanced numerical methods. This pertains to low numerical dissipation but also to low phase errors and low geometrical anisotropy errors because the flows of interest can quickly become geometrically complex. Minimizing numerical anisotropy (the imprinting of the grid) helps to preserve geometrical detail in these simulations. This can be accomplished by the use of approximately isotropic discretization methods [see the wave propagation method [START_REF] Leveque | Wave propagation algorithms for multidimensional hyperbolic systems[END_REF]], or simply by utilizing high-order methods. For example, Figure 4 shows a comparison of a spherical bubble undergoing radial collapse computed with the same mesh and with low-order (second-order HLLC) and high-order (fifth-order, weighted essentially nonoscillatory) numerical methods. An initially spherical bubble residing in a Cartesian cubic mesh implodes due to higher pressure outside the bubble; the results shown correspond to the time when the bubble reached its minimum size (see Tiwari et al. 2013 for details). Evidently, the high-order method preserves the sphericity of the bubble much better.

SOME RECENT RESULTS

In this section, we briefly discuss some recent applications of DIMs to various physical situations. Undoubtedly, the problem that has received most attention is shock-bubble interaction, followed by single bubble and cloud cavitation. Shock-bubble interaction has been studied for many decades, first in two dimensions and only more recently using 3D simulations. The substantial literature concerned with simulations of 2D shock-bubble interactions is not reviewed here due to space limitations, but readers are referred to [START_REF] Johnsen | Numerical simulations of non-spherical bubble collapse[END_REF] for some further references. More recently, [START_REF] Hejazialhosseini | Vortex dynamics in 3D shock-bubble interaction[END_REF] simulated a 3D shock-bubble interaction using the Kapila model. [START_REF] Meng | Numerical simulation of droplet aerobreakup[END_REF] also simulated a 3D shock-bubble interaction using a two-phase Euler model in the stripping breakup regime without what we call relaxation terms, but using the THINC interface sharpening method of [START_REF] Shyue | An Eulerian interface sharpening algorithm for compressible two-phase flow: the algebraic THINC approach[END_REF]. Figure 5 shows snapshots of these bubbles after the interaction with the shock. One can appreciate the complex breakup process and the formation of small-scale secondary bubbles, as well as the generation of vorticity.

Another situation of interest is the collapse of a bubble cloud near a wall (due to initial overpressure inside the bubbles). We highlight the work of [START_REF] Rossinelli | 11 PFLOP/s simulation of cloud cavitation collapse[END_REF] that has taken the Kapila model to new extremes by simulating the collapse of 15,000 bubbles near a wall. Figure 6 shows a snapshot of that simulation, which was possible thanks to algorithmic improvements recognized by the Association for Computing Machinery's 2013 Gordon Bell Prize. Finally, [START_REF] Tiwari | Growth-and-collapse dynamics of small bubble clusters near a wall[END_REF] performed simulations of a 50-bubble-cluster collapse using the Kapila model and Figure 6 A snapshot from the simulation of a 15,000-bubble-cluster collapse near a wall. Reproduced with permission from [START_REF] Rossinelli | 11 PFLOP/s simulation of cloud cavitation collapse[END_REF]. the interface compression of Section 3.5.1. The focus of this investigation was to understand the induced overpressure at the wall, where one expects the most damage.

Figure 7 shows the atomization of an oxygen liquid jet with a coaxial hydrogen gas flow in the absence of gravity and surface tension. This is a representative application of cryogenic fuel injection in the combustion chamber of space launchers. The figure compares the effect of phase change on the liquid oxygen dynamics. One observes that modeling phase change effects at the interphase of the oxygen and hydrogen fluids results in, given the same nominal conditions, very distinct spreading, concentration evolution, and interactions. Figure 8 considers the same flow model in the presence of both gravity and surface tension effects. Because the Mach number is very low, implicit time integration and preconditioning are required (Le [START_REF] Martelot | Liquid and liquid-gas flows at all speeds[END_REF].

Finally, we conclude with a recent application of the Kapila model and its PNE variant to the simulation of hyperelastic materials [START_REF] Favrie | Solid-fluid diffuse interface model in cases of extreme deformations[END_REF]. Additional equations are obviously needed to express the solid phase stress tensor, whereas the plastic transition is considered following [START_REF] Ndanou | Multi-solid and multi-fluid diffuse interface model: applications to dynamic fracture and fragmentation[END_REF]. The Kapila model enables the dynamic creation of interfaces, as a 0 ms 50 ms 100 ms 200 ms 300 ms 400 ms

Simulation of boiling flow with a liquid initially at saturation and containing three nucleation sites. Constant and uniform heat flux is imposed at the bottom wall. Vapor mass fraction contours are shown at 0, 50, 100, 200, 300 and 400 ms. Adapted with permission from Le [START_REF] Martelot | Liquid and liquid-gas flows at all speeds[END_REF].

T I M E

Impact of an iron spherical projectile onto an aluminum plate at 1,000 m/s. In addition to shock waves and interface motion, cracks appear and propagate dynamically due to the capability of the flow model to deal with the dynamic appearance of the interface in conjunction with plastic transformation. Adapted with permission from [START_REF] Ndanou | Multi-solid and multi-fluid diffuse interface model: applications to dynamic fracture and fragmentation[END_REF].

consequence of pressure relaxation effects, which, when combined with plastic transition, allows the modeling of dynamic fragmentation and crack propagation. Figure 9 shows an example of a high-speed impact of a solid projectile on a solid plate.

SUMMARY POINTS

1. DIMs combine the benefits of shock-/interface-capturing numerical methods, often based on Riemann solvers, with specialized formulations of the multicomponent inviscid compressible equations.

2. The fluids are allowed to interpenetrate across a small length scale while maintaining correct interface conditions.

3. Contrary to single-phase flow, the resulting models are not conservative in nature, and this necessitates special numerical methods.

4. Generalization to problems with surface tension, phase change, boiling, and cavitation has been implemented. More complex physical situations are still in development.

FUTURE ISSUES

1. It would be desirable to determine which treatment of nonconservative terms in DIMs' governing equations is best and, in particular, to understand which methods lead to converged solutions, both in 1D and in multidimensional problems. For example, it is hard to impose a given system of shock relations at the discrete level. This challenging problem is compounded by the potential nonlinearities present in the EOSs, which may complicate the stability of existing methods.

2. A comparative analysis should be performed to determine which sharpening/compression method produces the best results under a wide range of conditions or, for lack of a general answer, to determine which methods work best in different circumstances.

3. Robust solvers should be designed for flows at extremely high and low densities to ensure realizability of all fields in the presence of strong shocks and near vacuum conditions (cavitation).

4. Methods for low-Mach number flows should be improved. This requires the development of efficient preconditioning techniques for two-phase, DNS-like flow simulations.

5. Simple and well-posed (convex) EOSs should be developed in the presence of phase change. Because DIMs are hyperbolic, cubic EOSs have to be reconsidered, modified, and simplified in order to be useful in relaxation solvers (to ensure mechanical and thermodynamical equilibrium).

6. The existing methods should be extended to more complex fluids, e.g., non-Newtonian materials, elastic-plastic fluids with damage.

7. There is a growing need to develop well-defined test problems that possess unique solutions in order to validate the numerical methods. It is currently not routine in two-phase flow simulations to carry out convergence studies with grid refinement in multidimensional problems.

8. The multiscale nature of two-phase flows necessitates advanced and adaptive numerical methods: dynamic grid refinement and method order accuracy adaptation, efficient load balancing on multiprocessor supercomputers, and the development of hybrid methods that can blend averaged and resolved-interface formulations (e.g., for secondary atomization).

Figure 2

 2 Figure 2Simulation results without (a-c) and with (d-f ) compression. (a,d ) Volume fraction comparison of shock-bubble interaction, adapted with permission from[START_REF] Shukla | An interface capturing method for the simulation of multi-phase compressible flows[END_REF]. (b,e) Schlieren of density for shock interaction simulation with a R22 gas bubble, adapted with permission from[START_REF] Shyue | An Eulerian interface sharpening algorithm for compressible two-phase flow: the algebraic THINC approach[END_REF]. (c,f ) Unstructured grid simulation results of shock interaction with a gas bubble. In each row, the same resolution was used in corresponding left and right simulations, with the shock traveling from left to right, but unfortunately the conditions and resolution varied from the top row to the bottom row.

  Figure 3

  Figure 5Shock-bubble interaction simulations from (a)[START_REF] Hejazialhosseini | Vortex dynamics in 3D shock-bubble interaction[END_REF], with the shock going from left to right, and (b) Meng (2016), with the shock going from right to left and with isosurfaces of positive (red ) and negative (blue) azimuthal vorticity indicated. Adapted with permission from[START_REF] Hejazialhosseini | Vortex dynamics in 3D shock-bubble interaction[END_REF] and[START_REF] Meng | Numerical simulation of droplet aerobreakup[END_REF], respectively.

Figure 7

 7 Figure 7 Simulation of liquid oxygen fragmentation (Y 1 ) by a coaxial hydrogen gaseous co-flow (Y 2 ) in conditions of cryotechnic engine combustion restart. Mass fraction contours shown (a) without phase change Y 1 , (b) with phase change Y 1 , (c) with phase change Y 1 + Y 2 , and (d ) with phase change Y 2 . Readers are referred to Chiapolino et al. (2017a) for details.

  discussed efficient techniques for solving the volume fraction equation with compression for large U o . The second method is also an Eulerian-sharpening algorithm referred to as the tangent of hyperbola for interface capturing (THINC) (seeShyue &

	3.5.2. Interface sharpening technique.

ACKNOWLEDGMENTS

We thank Professors K.-M. Shyue, T. Colonius, E. Johnsen, P. Koumoutsakos, and N. Adams, as well as Dr. J.C. Meng, for sharing some of their figures.

DISCLOSURE STATEMENT

The authors are not aware of any biases that might be perceived as affecting the objectivity of this review.