CWENO Finite-Volume Interface Capturing Schemes for Multicomponent Flows Using Unstructured Meshes

In this paper we extend the application of unstructured high-order finite-volume central-weighted essentially non-oscillatory (CWENO) schemes to multicomponent flows using the interface capturing paradigm. The developed method achieves high-order accurate solution in smooth regions, while providing oscillation free solutions at discontinuous regions. The schemes are inherently compact in the sense that the central stencils employed are as compact as possible, and that the directional stencils are reduced in size, therefore simplifying their implementation. Several parameters that influence the performance of the schemes are investigated, such as reconstruction variables and their reconstruction order. The performance of the schemes is assessed under a series of stringent test problems consisting of various combinations of gases and liquids, and compared against analytical solutions, computational and experimental results available in the literature. The results obtained demonstrate the robustness of the new schemes for several applications, as well as their limitations within the present interface-capturing implementation.


Introduction
Multi-component multi-phase compressible flows are encountered in several scientific fields and due to the complexity of the interaction of different states of matter a detailed understanding of them is required for many applications. The choice of schemes to simulate these interactions are far neither unique nor obvious, due to several factors such as the level of fidelity, the computational resources, and due to the mechanism of the interface that separates two fluids (e.g. liquid and gas), since the thermodynamic properties of the interface might be different from the two separate states, that can obey different equation of state (EOS). Most of the schemes in an Eulerian framework belong to the class of sharp inter-face type methods (SIM) or the class of diffuse interface methods (DIM). The SIM class of methods treat the material interface between the fluids as infinitely thin, or sharp and several successful methods have been developed including the volume-of-fluid [28], the front-tracking [55], and the level-set methods [40]. The DIM class of methods is treating the interface as a non-zero thickness based on the idea originally formulated by Rayleigh [45] and van der Waals [69] who developed gradient theories for predicting the thickness of the interfaces based on thermodynamic principles. The DIM allow artificially the two states to develop some mixing, and therefore a thermodynamic state of this mixture is required as the interface evolves. Among the common DIM models are the five-equation model of Kapila's et al. [35], which is tailored for two-phase immiscible non-reacting flows which was also adopted by Allaire's et al. [1]. The complete Baer-Nunziato's [6] seven equation model that includes non-equilibrium effects and the unified hyperbolic formulation of Godunov-Peshkov-Romenski (GPR model) [24,42] that can treat traditional fluid and solid problems within one framework which is the ultimate goal of several highly sophisticated and elegant advanced diffuse interface models [8,14,22,23,29,38,70,71].
In the present study the five-equation model of Kapila et al. [35] is employed in the form adopted by Allaire et al. [1] where an isobaric closure law is deployed that can simulate two fluids with arbitrary EOSs. This form has been extended by numerous researchers to include viscosity, capillary effects, and has been also extended to unstructured meshes by Chiapolino et al. [11], Price et al. [43], Faucher et al. [21] and Cheng et al. [10]. The simplified and more compact five-equation model of Allaire et al. has found many applications as listed in [10,11,13,21,32,36,43,47,72,73].
The hyperbolic character of the PDEs involved in this context, the simplicity of this type of DIM and their ease of implementation are the primary reasons for the selection of this five-equation model. For engineering applications unstructured meshes offer significant advantages when dealing with optimisation of complicated geometries that require computational workflows with rapid and highly automated mesh generation. In order for the DIM class of methods to work correctly the spatial accuracy should be sufficient to resolve the interfaces correctly without an excessive diffusion across them and this requirement can be achieved by either increasing the spatial resolution at the interface regions or resorting to highresolution numerical methods. Additionally the numerical methods should be non-oscillatory since very strong gradients across material interfaces can result in spurious oscillations and a blown-up simulation. There are several high-order high-resolution numerical methods in the unstructured finite-volume framework such as the Weighted Essentially Non-Oscillatory (WENO) [17,63,75], Central Weighted Essentialy Non-Oscillatory (CWENO) [19,61], Multidimensional Optimal Order Detection (MOOD) [15,20], Monotone Upstream Scheme for Conservation laws (MUSCL) [52,57] and the Edge Based Reconstruction WENO (EBR-WENO) [7].
In this study however we will be focusing on the application of the CWENO type of schemes as developed in [19,61] since these methods were found to be more robust and significantly more computational efficient compared to the original WENO type of schemes for unstructured meshes, that have previously been successfully applied by Dumbser et al. [18] to the full seven equation Baer-Nunziato model. The improved robustness and efficiency is provided by the low-order polynomials associated with the directional stencils and their reduced size compared to the high-order polynomial associated with the central stencil. The reader is referred to [19,61] for an overview of the methods. All the schemes/models are developed in the open source community UCNS3D solver [68], and we assess their performance in terms of for a series of stringent 2D and 3D test problems. The paper is organized as follows. In Sect. 2 we introduce the formulation of the governing equations employed for this study, followed by the numerical framework used to describe the highorder finite-volume framework for unstructured meshes, the reconstruction process for the CWENO schemes, while describing the chosen fluxes and temporal discretisation employed in Sect. 3. In Sect. 4 the numerical results obtained for all the test problems are presented and compared against analytical, reference or experimental solution whenever possible. Finally, the last section describes the conclusions drawn from this study.

Governing Equations
The quasi-conservative five-equation model of Allaire et al. [1] is considered in this study for inviscid compressible multicomponent flows. For two fluids this involves two continuity equations Eqs. (1) and (2), a momentum equation per dimension Eq. (3), an energy equation Eq. (4), and the non-conservative advection equation of the volume fraction of one of the two fluids Eq. (5) as given below: ∂ρu ∂t where ρ is the density, u = (u, v, w) T is the velocity, p is the pressure, E is the total energy and a is the volume fraction. The widely use stiffened gas EOS is employed for closing the five-equation model. It has been primarily selected due to its application for flow problems involving gases, liquids and solids. The stiffened gas EOS characterises each fluid pressure as: where π ∞,i ≥ 0 is a reference pressure, and will be set to π ∞ = 0 for gases. The total mass and ρ being given by the following equations: where is the internal energy, with ρ = E − 1 2 ρuu. The EOS of the mixture reads Finally, the non-conservative volume fraction advection Eq. (5) is rewritten in a mathematically equivalent form as introduced by Johnsen and Colonius [32]: 3 Numerical Framework

Spatial Discretization
Consider a 3D domain Ω consisting of conforming tetrahedral, hexahedral, prism, and pyramid cells each one of them indexed by a unique mono-index i, and the governing equations of the five-equation model written in vector form as follows: where U = U(x, t) is the vector of conserved variables and the volume fraction of one species, F n is the non-linear flux in the direction normal to the cells interface as given below: where u n is the velocity normal to the bounded surface area, defined by u n = n x u+n y v+n z w. The source term s is with regards to the term a 1 ∇ · u of Eq. (5). Following the approach of Johnsen and Colonius [32] the source term is numerically approximated as surface integral, rather than a volume one, while using the same velocity estimate as the one used for the evaluation of the fluxes as shown below: Integrating Eq. (13) over the mesh element i using a high-order explicit finite-volume formulation the following equation is obtained that incorporates the source term as previously defined: where U i is the volume averaged vector of variables and F n ij is a numerical flux function in the direction normal to the cell interface between a considered cell i and one of its neighbouring cells j. N f is the number of faces per element, N qp is the number of quadrature points used for approximating the surface integrals, |S i j | is the surface area of the corresponding face, and U n i j,L (x i j,α , t) and U n i j,R (x i j,α , t) are the high-order approximations of the solutions for cell i and cell j respectively. α corresponds to different Gaussian integration points x α and weights ω α over each face. a n i,1 corresponds to the volume averaged volume fraction of cell i at time level n. The volume, surface and line integrals are numerically approximated by a suitable Gauss-Legendre quadrature.

Reconstruction
A high-order polynomial p i (x, y, z) of order r can be built that provides r + 1 order of accuracy for a cell i , by ensuring that it has the same average as a general quantity U i . This can be formulated as: The present polynomial reconstruction is based upon the approaches of [53,63], that have been applied to smooth and discontinuous flow problems [2][3][4][5]20,46,49,50,57,[59][60][61][62][64][65][66][67], and only the key-components will be presented herein and the reader is referred to [53,63] for further details. For the reconstruction, a transformation from physical space to a reference space as introduced by Dumbser et al. [16,17] is employed to reduce the scaling effects that occur at unstructured meshes consisting of elements of different shape and size. In particular we employ the decomposition strategy defined in [63,64], where each element is decomposed into triangular (2D) or tetrahedral (3D) elements and using one of the decomposed elements as the reference element for transforming to the new system of coordinates. Let vx i j , j = 1, 2, . . . J i be the vertices of the considered 3D general element. Non tetrahedral elements are decomposed into tetrahedrals and one of them is chosen with w 1 = (x 1 , y 1 , z 1 ), w 2 = (x 2 , y 2 , z 2 ), w 3 = (x 3 , y 3 , z 3 ) and w 4 = (x 4 , y 4 , z 4 ) being its four vertices. The transformation from the Cartesian coordinates x, y, z into a reference space ξ, η, ζ is given by the following equations: with the Jacobian matrix given by: Using an inverse mapping the element V i can be transformed to the element V i in the reference co-ordinate system as: while the spatial average of the conserved variable U i does not change during transformation The reconstruction polynomial uses a central stencil S 1 , which is built by recursively adding neighbouring elements, consisting of M + 1 cells including the considered cell i. We employ the stencil based compact algorithm (SBC) introduced in [58], which is characteristed by a low computational cost and increased robustness in regions of poor quality meshing. In the present study we employ M = 2K , for enhanced robustness as reported in several previous studies [15,30,58], where K is the total number of polynomial coefficients given by: where d ∈ [2,3] is the space dimensions. The central stencil S 1 is given by where the index m refers to the local numbering of the elements in the stencil with the element with index 0 being the considered cell i, and the index c referring to the stencil number (in case of multiple stencils) where for the central stencil c = 1. All the cells of each stencil are transformed in the reference space S c i , where the r th order reconstruction polynomial is an expansion over local polynomial basis functions φ k (ξ, η, ζ ) given by: where U 0 corresponds to the vector of conserved variables at the considered cell i, and a k are the degrees of freedom of the polynomial. The degrees of freedom a k for the polynomial for each cell m are obtained by ensuring that the cell average of the reconstruction polynomial p(ξ, η, ζ ) is equal to the cell average of the solution U m : The conservation condition in Eq. (18) is such that requires the basis functions to have a zero mean value over the considered transformed cell V 0 . This can be achieved by using hierarchical orthogonal basis functions defined on a unit element in reference space. For triangular and tetrahedral elements hierarchical orthogonal basis functions can achieve this requirement. However for arbitrary shaped quadrilaterals, pyramids, prisms and hexahedrals they do not transform to a unit-element exactly, and therefore their basis functions need to be constructed in a way that the conservation condition in Eq. (18) is satisfied. The basis functions employed φ k for all the elements in the stencil are defined as follows and satisfy this requirement: and in the present study ψ k are Legendre polynomials basis functions. Denoting the integrals of the basis function k over the cell m in the stencil, and the vector of right-hand side by A mk and b respectively as given by the equations for degrees of freedom a k can be rewritten in a matrix form as: The resulting linear system is solved by a QR decomposition based on Householder transformation while using a Moore-Penrose pseudo-inverse of A mk which is only computed once at the beginning of the simulation as detailed in [58].

CWENO Scheme
The CWENO scheme employed follows the implementation of Tsoutsanis and Dumbser [61] and is the combination of an optimal (high-order) polynomial p opt with lower-order polynomials. The reconstruction for the optimal (high-order) polynomial uses an expanded central stencil, (suitable for the desired polynomial order) while the reconstruction for the lower-order polynomials employs the compact directional stencils. When the variation of the solution is smooth across all the stencils, the optimal polynomial is recovered and therefore the desired-order of accuracy is obtained. At the presence of discontinuous data at least one of the directional stencils associated with the lower-order polynomials could contain smooth data, hence it is going to have the largest influence in the reconstruction. All the polynomials involved satisfy the requirement of matching the cell averages of the solution, therefore they are solved with the same constrained least-squares technique. The directional stencils employ the Type3 definition which includes one directional stencil per element face as detailed in the work by Tsoutsanis [58]. The optimal polynomial is defined as follows: where s is the stencil index, with s = 1 being the central, s = 2, 3, being the directional, s t being the total number of stencils, and λ s being the linear coefficients for each stencil, whose sum is equal to 1. The p 1 polynomial is not computed directly, but computed by subtracting the lower-order polynomials from the optimum polynomial as follows: The CWENO reconstruction polynomial is given as a non-linear combination of all the polynomials in the following manner: where ω s correspond to the non-linear weights assigned to each polynomial, and in regions with smooth data ω s ≈ λ s , hence obtaining the high-order approximation from the central stencil, and in regions of discontinuous solutions the reconstructed solution will be mostly influenced from the lower-order polynomials of the directional stencils. whereã k are the reconstructed degrees of freedom; and the non-linear weight ω s is defined as: The smoothness indicator SI m is given by: where β is a multi-index, r is the polynomial's order, λ m is the linear weight. The value set to prevent division by zero of = 10 −6 is used, with b = 4 and D being the derivative operator. The smoothness indicator is a quadratic function of the degrees of freedom (a s k ) and Eq. (33) can be rewritten as: where the oscillation indication matrix OI kq is given by: which can be precomputed and stored at the beginning of the simulation. For the directional stencils and their corresponding polynomials we employ r = 1 to obtain 2nd-order of accuracy, and any arbitrary order of accuracy for the polynomial associated with the central stencil.
The linear weights are computed by firstly assigning the non-normalised linear weight for the central stencil λ 1 an arbitrary value, and then normalising this as follows: with the linear weights associated with lower-order polynomials being equal and provided by the following expression: where s t is the total number of stencils. From this point forward the order of the scheme will be defined by a number next to the type of the scheme, such as CWENO3 and WENO4 corresponding to a 3rd-order CWENO scheme and a 4th-order WENO scheme respectively.

Fluxes Approximation and Temporal Discretisation
For the inviscid fluxes the approximate Harten-Lax-van Leer-Contact (HLLC) Riemann solver of Toro [54] is employed, unless otherwise stated. The temporal discretisation employs the 3rd-order explicit Strong Stability Preserving (SSP) Runge-Kutta method [26] which is stable for C F L ≤ 1. All the volume/surface/line integrals are approximated by Gaussian quadrature rule suitable for the order of polynomial employed. It has to be noted that the reconstruction in the present study is carried out with respect to the conserved variables or primitive variables by transforming to either one of them only during the reconstruction process, while always transforming to conservative variables prior to the fluxes calculation. All the schemes developed are implemented in the open-source UCNS3D CFD code [68] using object-oriented Fortran 2003, and employing MPI message passing interface (MPI), and the Open Multi-Processing (OpenMP) application programming interface (API). The reader is referred to [56,67] for more details on implementation and performance benchmarks.

Applications
We present the numerical simulations employed to assess the performance of the CWENO interface capturing schemes for the solution of the inviscid compressible Euler equations for multicomponent flows. Several benchmark test problems have been performed.

Multi-Species Convergence Test
For verifying the designed order of spatial accuracy for the numerical schemes developed, a multi-species advection test similar to the one employed by Wong and Lele [73] is used. In this test a smooth volume fraction initial profile of two gases is advected for one period in a periodic computational domain. The initial condition is given by: The 2D computational domain [0, 1] 2 consists of arbitrary unstructured triangular elements of 10, 20, 40 and 80 edges per side resolution as shown in Fig. 1, and the simulation is run for a time of t f = 1. The two gases selected are nitrogen and helium with specific heats 1.4 and 1.66 respectively. The numerical errors e L 2 and the e L ∞ are computed as follows: where U c x, t f and U e x, t f are the computed and exact solutions at the end of the simulation t f = 1.0. The exact solution U e x, t f being given by the initial condition itself at t = 0. The simulations were performed using C F L = 0.1 in order to ensure a sufficiently small time-step size so that the schemes are not restricted by the time-discretisation for achieving the designed order of spatial accuracy. The wallclock time per simulation is normalised with respect to the fastest time taken for a simulation, which in this setup is provided by the CWENO3 scheme on the coarsest mesh. From the obtained convergence orders of the schemes as shown in Table 1 it is clear that the schemes achieve their designed order of convergence for both e L ∞ and e L 2 . The desirable feature of the CWENO variant is that it comes with significant lower computational costdue to the reduction of the size of the directional stencils-compared to the WENO schemes. For this test the CWENO schemes require approximately 60-80% of the time taken by a WENO scheme of the same spatial order, a figure which is expected to improve when larger mesh sizes, and 3D test problems are deployed. It can be noticed that all the schemes achieve their designed order of accuracy, with the CWENO schemes being requiring approximately (60-80%) of the time taken by an equivalent WENO scheme

Isolated Material Interface
For assessing the non-oscillatory properties of the considered schemes, the advection of a sharp material interface within a periodic domain is considered. A sharp material interface is frequently encountered in several multicomponent flows, hence the robustness of the proposed methods is of paramount importance, for the successful deployment of these schemes in multicomponent flows. The material interface is advected with constant velocity across the computational domain, and the pressure and temperature is also constant across the interface. The initial conditions are given by: The 2D computational domain [0, 1] 2 consists of a mixed-element mesh as shown in Fig. 1 consisting of triangular and quadrilateral mesh elements of 80 edges per side resolution, and the simulation is run for a time of t f = 2. Firstly a comparison between 3rd-order CWENO3 and WENO3 is made to assess the non-oscillatory properties of the schemes and as it can be noticed from the obtained results in Fig. 2 the CWENO3 is significantly less oscillatory near the material interface. The relative coarse resolution of the present mesh, is intentionally selected to highlight one of the key differences between the two schemes, which is the compactness of their directional stencils. In the CWENO the directional stencils correspond to lower-order polynomials compared to WENO schemes and therefore are significantly more compact, therefore resulting in improved robustness in such situations. Secondly the influence of two types of reconstructions are explored with respect to the conserved and primitive variables. As it can be seen in Fig. 3, the solutions obtained with the primitive variable reconstruction are free from any oscillations in the pressure and velocity solution. It is well documented [31][32][33][34] that primitive variable reconstruction is needed to prevent spurious oscillations in problems where γ is not constant, since reconstruction in conservative or characteristic variables suffer from oscillations across them. Finally the CWENO scheme with primitive variables employed provide the correct solution with no oscillations, and as seen in Fig. 2 Plots of density a 1 ρ 1 for the isolated sharp-material interface at t = 2 obtained with WENO3 and CWENO3 schemes, and compared with the reference exact solution. It can be noticed that the the WENO3 is producing some oscillations near the material interface, while they are absent from the CWENO3 obtained solution Fig. 3 Plots of density a 1 ρ 1 , pressure and velocity for the isolated sharp-material interface at t = 2 obtained with the CWENO3 scheme with primitive and conservative variables reconstruction, and comparison with the exact solution. The primitive variable reconstruction is producing an oscillation-free solution for both pressure and velocity Fig. 4 Plot of density a 1 ρ 1 , pressure and velocity for the isolated sharp-material interface at t = 2 obtained with several CWENO schemes where it can be noticed that all the schemes provide the correct solution with no oscillations, and normalised pressure and velocity error at the t = 2 where it can be seen that the minute oscillations are close to machine precision Fig. 4 the minute oscillations in normalised pressure and velocity errors are close to machine precision. Therefore unless otherwise stated from this point onwards the reconstruction with respect to primitive variables is going to be employed.

Gas-Liquid Riemann Problem
The 2D equivalent of the well established gas liquid Riemann problem of Cocchi et al. [12] is deployed in this study. It is an ideal test problem for exposing the performance of the developed numerical schemes in for stiffened gas EOS. The left state of the problem consists of highly compressed air, and the right state is water at atmospheric pressure. The initial conditions for this problem following the non-dimensionalisation of Cocchi et al. [12] are: 241, 0, 0, 2.753, 1.4, 1) , The specific heats for air and water are 1.4 and 5.5 respectively, with the π ∞ = 1.505 for water. The 2D computational domain [−1, 1] 2 is discretised by a mixed-element unstructured mesh similar to the one shown in Fig. 1 which consists of arbitrary triangular and quadrilateral elements of 100 edges per side resolution, and the simulation is run for a time of t f = 0.2. Two schemes are employed for this test problem a WENO3 and a CWENO5. The CWENO5 scheme employs a higher-order polynomial for the central stencil compared to the WENO3 order scheme, and a lower-order polynomial (and stencil size) than the WENO3 order scheme. As it can be noticed from the obtained results shown in Fig. 5 while both schemes provide results in good agreement with the exact solution, the CWENO5 order scheme is less oscillatory in this coarse grid resolution due to the compactness offered by the smaller directional stencils. Therefore from this point onwards only the CWENO schemes will be employed in the following test problems.

Gas Bubble in Water
A problem that has been frequently used [25,39] to assess the performance of various numerical methods is considered. This problem involves very high-pressure due to the collapse of a gas bubble in a liquid. These types of problems are typically found in applications such as fuel injectors, naval propulsion systems and shockwave lithotripsy. In this study a strong shock M sh = 1.72 is travelling through water and moving towards an air bubble. The initial condition is given by: The evolution of the bubble dynamics obtained with the present schemes is illustrated in Fig. 6 and in Fig. 7 where all the expected features including the reflected rarefaction wave, the water jet, the blast wave and the secondary jets are well captured and are similar to Plots of density, pressure, velocity and volume fraction for the gas-liquid Riemann problem at t = 0.2 obtained with WENO3, and CWENO5 schemes and compared with the exact solution. It can be noticed that the CWENO5 order scheme is less prone to oscillations compared to the WENO3 order scheme the ones obtained by others Nourgaliev et al. [39] and Goncalves et al. [25]. The maximum pressure obtained in the computational domain is 7.7G Pa, which is inline with the findings of other studies where higher pressures are seen in fine grid resolutions with Nourgaliev et al. [39] reporting 10.1G Pa at a grid resolution of e c ≈ D b /800, and lower peak pressures at lower grid resolutions with Goncalves et al. [25] reporting a peak pressure of 7.8G Pa at a grid resolution of e c ≈ D b /300. Similarly the peak temperature at the time of lowest volume obtained is ≈ 15, 670K which is within the range obtained from the other studies [25,39].
The time evolution of the pressure along the centreline of the computational domain can be seen in Fig. 8, where it can be noticed that the first peak is associated with the blast wave at t ≈ 3.9 µs and the second peak a later time t ≈ 4.9 µs occurs when the blast wave collides with the bubble fragments as also reported by Goncalves et al. [25]. Finally the time history of the non-dimensional volume of the gas bubble is illustrated in Fig. 8, where it can be noticed that the results obtained are in good agreement with the results of Nourgaliev et al. Fig. 6 Plots of the Mach number (contours) and volume fraction (lines) for the gas bubble collapse test problem at different instants. It can be noticed that all the expected features are captured using the present CWENO5 scheme . 7 Plots of the density contours (top), volume fraction (middle), and density gradient magnitude (bottom) for the gas bubble collapse test problem at different instants [39] capturing the correct compression slope, and rebound at late times. However due to the coarser grid-resolution employed in this study the maximum compression of the bubble is lower compared to the one obtained by Nourgaliev et al. [39].

Water Column in Air
The shock wave interaction with a cylindrical water column test case of Xiang and Wang [74] is employed in this study, and in essence it is the opposite of the gas-bubble collapse in water problem seen previously. In this test a water droplet is surrounded by air, where a shockwave at M sh = 2.4 is moving towards the water bubble. The subject test serves as an ideal computational test-problem for development and assessment of numerical methods for interface dynamics of compressible multi-fluid flows with several practical applications including the liquid jet atomization in a scramjet engines, combustion, and supernova explosion. We are considering two variants of the water column in this test one without an air cavity, and one with an air cavity as performed by Xiang and Wang [74]. The setup of the test problem can be seen in Fig. 9, where for the water column without the air cavity r = 0, while with air cavity r = 3.6 mm. The initial condition is given by: 3.85, 567.3, 0, 0.664 · 10 6 , 1.4, 0, 0 , for Post-shock 1.2, 0, 0, 0.101 · 10 6 , 1.4, 0, 0 , for Pre-shock 1000, 0, 0, 10 5 , 6.12, 0.343 · 10 9 , 1 , for Water (44) The domain is x ∈ [0, 0.2208], y ∈ [0, 0.1152], the water column is centred at (0.0576 m, 0.0576 m), and the interface between shocked and unshocked regions placed at x = 0.05 m. Non-reflecting boundary conditions are prescribed at the left and right boundaries of the domain, while an Euler slip-wall is prescribed for the top and bottom boundaries. The domain is discretised by an quadrilateral mesh with approximately 0.89 million cells and an average edge length of e c ≈ D b /192, where D b is the diameter of the water column/or the outer diameter when a cavity is included D b = 96 mm. The CWENO5 scheme was employed for this test problem, and the simulation is run until the non dimensional t * = tu/D b = 10 where u corresponds to the initial velocity behind the shock. The obtained computed density gradient magnitude for both variants are in good agreement with the numerical results of Xiang and Wang [74] as shown in Fig. 10. All the key flow features (reflected expansion wave, the transmitted wave, the Mach stem, and re-circulation zone) are correctly captured Fig. 9 Schematic diagram of the setup for the shock wave cylindrical water column test problem (left), and the corresponding mesh refinement used for the interaction zone (right) Fig. 10 Contour plots of density gradient magnitude of the computed solution of the shock wave interaction with a water column with a cavity (top) and without a cavity (bottom), at various instants t * = 0, t * = 0.8, t * = 1.62, t * = 3.02 (left to right). It can be noticed that all the expected flow features are captured by the present CWENO scheme as shown in comparison with the experimental results for the water column without cavity of Sembian et al. [48] in Fig. 11.
As expected at late times the flow separates for both variants as seen in Fig. 12, where the water column without the cavity is compressed continuously while it gets flattened downstream as documented by Xiang and Wang [74] and Meng and Colonius [37]. For the water column with cavity the transverse jet formation and the associated Richtmyer-Meshkov instability occurring when the water is driven into the air cavity are well resolved, and in agreement with the study of Xiang and Wang [74].
Finally, the pressure along the centreline of the water column without the cavity at time t * = 0.8 is plotted against the computational results of Xiang and Wang [74] in Fig. 13 and a good agreement is achieved. It has to be stressed that the expected negative pressure drop is not as pronounced as the study of Xiang and Wang [74] since their grid resolution was finer with an edge length of of e c ≈ D b /333. However, the negative pressure is correctly captured due to the reflected expansion wave, which is attributed to the cohesive forces that hold liquids together and can withstand negative pressure [9,41]. For the water column with the air cavity, a time-history of the maximum pressure present in the computational domain is plotted against the computational results of Xiang and Wang [74] in Fig. 13, where it is evident that a good agreement is achieved although the peak pressure is slightly lower due to the lower grid resolution as expected. However the timing of the maximum peak pressure agrees with the results of [74], where the maximum peak pressure is seen at the time of impact of the transverse jet with the inner cavity wall downstream.

Helium Bubble Shock Wave
The interaction of a weak shockwave in air and a helium bubble is considered in 2D and 3D. Several variations of this test problem have been widely used [13,32,72] for assessing the performance of several techniques for multicomponent flow modelling, and is based on the experimental setup by Haas and Sturtevant [27]. A bubble of diameter D b = 5 cm, is placed within an air filled shock tube. The bubble consists of helium and air of 28% mass concentration. A shockwave moving from right to left as shown in Fig. 14 of the setup impacts the bubble contaminated by the surrounding air. The specific heats of 1.4 and 1.66 are used for air and helium, respectively, and the initial condition is given by: (a 1 ρ 1 , a 2 ρ 2 , u, v, p, a 1 The computation domain is discretised by a mixed-element unstructured mesh consisting of quadrilateral and triangular elements in 2D, and arbitrary hexahedrals and prism in 3D with  Plot of the definition of the interfaces and of the evolution of the position of these interfaces using a CWENO5 scheme on the finest 2D mesh and 3D mesh. The evolution of the interfaces position is in good agreement with the results of Terashima and Tryggvason [51] and Quirk and Karni [44] the shock-bubble interaction region being refined as shown in Fig. 14. A coarse, medium, and fine mesh are employed for 2D with an average element edge length in the shock bubble interaction zone of e c ≈ D b /35, e m ≈ D b /120, e f ≈ D b /400 respectively, and an edge length of e c ≈ D b /50 for the 3D mesh. A slip-wall boundary condition is used at the top and bottom boundaries of the domain, while inflow and outflow boundary conditions are prescribed at the right and left of the domain respectively. A CWENO5 scheme is employed and the simulation is run until t = 1000 µs. From the results obtained as shown in Fig. 16, it can be noticed that the time evolution of the bubble is correctly captured including the formation of a jet and a vortex ring at late times, and is in agreement with the results obtained by [13,32,72] qualitatively. Due to the high-order CWENO scheme employed, the Kelvin-Helmholtz instabilities that develop at the interface of the helium bubble are more pronounced as the 2D grid resolution is increased as shown in Fig. 17, whereas for the 3D setup the resolution employed is not sufficient to capture these instabilities. as shown in as shown in Fig. 18. The symmetry of the computed solution is lost at late times for the 2D medium and fine meshes, indicative of multi-dimensional reconstruction nature of the framework employed and of the arbitrary unstructured elements in the refined shock-bubble interaction zone. As the bubble evolves, three distinct interfaces identified are the jet ( ji), the upstream (ui) and the downstream (di) as shown in Fig. 15. From the space-time diagram of these distinct interface positions, the predicted locations are in good agreement with the reference results of Terashima and Tryggvason [51] and Quirk and Karni [44]. It needs to be highlighted that the results have been non-dimensionalised with respect to the diameter of the bubble at the time that the shock hits the bubble [44,51,72].  The results obtained in terms of the averaged velocities of the jet, upstream and downstream interface are in good agreement with the experimental results of Haas and Sturtevant [27], and the computational results of Coralic and Colonius [13] as shown in Table 2. Due to the resolution employed for the 3D setup of the problem, a slower merging of the jet and the downstream interface is seen compared to the 2D setup.

Conclusions
This paper extends to applicability of CWENO schemes on unstructured meshes to compressible multi-material flows. The schemes manage to achieve high-order of accuracy, resolve the material interfaces and finer structures while maintaining their essentially non-oscillatory character. Switching to primitive variable reconstruction was needed to remove the oscillations that appear across material interfaces. A series of stringent two-and three-dimensional test problems were used to verify the accuracy, robustness and computational efficiency of the schemes, and compared with analytical, experimental and computational results. Future development will concern the expansion of the CWENO schemes to the complete seven equation model of Baer-Nunziato's [6] and the unified hyperbolic formulation of Godunov-Peshkov-Romenski (GPR model) [24,42]. appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.