Applied Mathematics and Computation

In this paper the relaxed, high-order, Multidimensional Optimal Order Detection (MOOD) framework is extended to the simulation of compressible multicomponent ﬂows on un-structured meshes. The diffuse interface methods (DIM) paradigm is used that employs a ﬁve-equation model. The implementation is performed in the open-source high-order unstructured compressible ﬂow solver UCNS3D . The high-order CWENO spatial discretisation is selected due to its reduced computational footprint and improved non-oscillatory behaviour compared to the original WENO variant. Fortifying the CWENO method with the relaxed MOOD technique has been necessary to further improve the robustness of the CWENO method. A series of challenging 2-D and 3-D compressible multicomponent ﬂow problems have been investigated, such as the interaction of a shock with a helium bubble, and a water droplet, and the shock-induced collapse of 2-D and 3-D bubbles arrays. Such problems are generally very stiff due to the strong gradients present, and it has been possible to tackle them using the extended MOOD-CWENO numerical framework .


Introduction
In most engineering and scientific applications involve compressible flows which may exhibit a variety of flow components, materials, and species, giving rise to very complex phenomena that need to be adequately understood for design of products in the industry, and experiments in academia.Examples in which a wealth of compressible flow features is displayed include advanced aerospace applications like hypersonic scramjet engines, where the interaction of shock waves with fuel droplets occurs in a highly-compressible and turbulent regime (shock-driven multiphase instabilities) [1] and they are ubiquitous in laboratory plasma physics and astrophysics, with prime examples being inertial confinement fusion and supernova explosions (Rayleigh-Taylor and Richtmyer-Meshkov instabilities) [2] .
This work attempts to add more knowledge in this field by formulating a series of numerical tests using the Eulerian code UCNS3D .In contrast to the more standard treatment in Eulerian codes, this work is exclusively using what is called diffuse interfaces , a concept first put forth by Rayleigh [3] and van der Waals [4] .Typically, fluid interfaces in a hydrodynamics code are infinitesimally thin.This technique corresponds to the type of linear stability analysis in which the various flow features can be decomposed into a series of discrete Fourier modes.Our work, making use of the diffuse interface method (or DIM for short), does away with this approximation given its inherent limitation.Sharp interfaces are adequate for capturing large perturbational modes, but fail to resolve the finer modes (i.e.turbulent structures) developing during mixing across interfaces.Hence, this approximation neglects the onset of minute, non-linear effects.A diffuse interface has a finite thickness, and can therefore capture smaller perturbational modes that manifest along the interface of two mixing components of the flow field.One may consider the so-called Wentzel-Kramers-Brillouin, or WKB approximation, developed in 1926 [5][6][7] as an analytical equivalent.By using the DIM in this study, it thus becomes possible to capture more accurately smaller scales in the inertial range of the turbulent cascade, and consistently with the respective thermodynamic states at each side of the interface, which fundamentally improves the predictive capabilities of our solver to capture complex effects, while it is programmatically simple to implement.
The problems we investigated in this study are all multicomponent, and given the complexity of each case, the equations of state (EOS) employed may vary arbitrarily.For this reason, our solver, UCNS3D , has implemented a fiveequation model first proposed by Allaire's et al. [8] ; this model assumes an isobaric-closure which enables the code to use of different EOS per component.Therefore, the model offers certain advantages, as outlined in previous work [9,10] .
Nevertheless, it is possible, and desirable to extend the capabilities of UCNS3D with more advanced and complete models in the future, and our insights from this work will pave the way for future code development.The reader might take interest in the work of Maltsev et al [11] for a comprehensive review of more complicated models.Allaire's et al. model is compact and relatively uninvolved, finding good use in a multitude of applications as it can be seen for instance in Schmidmayer et al. [12] , Kokkinakis et al. [13] , Coralic and Colonius [14] , Wang et al. [15] , Wong and Lele [16] , Johnsen and Colonius [17] , Chiapolino et al. [18] , Price et al. [19] , Cheng et al. [20] , Faucher et al. [21] and references therein.The model by Allaire et al. [8] has a relatively straightforward implementation, a key factor considering the complexities involved in an unstructured framework.Additionally, this model is well suited for wave equations which are hyperbolic.
While diffuse interfaces are greatly expanding the code's potential to capture finer turbulent structures, the have the propensity to unphysically diffuse over time.Several measures are required to control this diffusion.Such countermeasures involve complementing the DIM with interface sharpening techniques or adaptive mesh refinement, as noted in Maltsev et al. [11] .Alternatively, improvements in the resolution sharpness of the interfaces present in the computational domain as well as suppression of any excessive diffusion across them are possible by employing high-order, spatially high-resolution, oscillation-free numerical methods.Nonetheless, one of the challenges that code development is facing when employing high-order numerical methods using the DIM approximation is the presence of very strong gradients of the flow variables in interfaces across which state of matter is different.Therefore some additional enhancement of the selected numerical methods might be necessary.
Most state-of-art highly accurate numerical schemes use some form of adaptive criteria to adjust the solution's spatial accuracy when sharp gradients are predicted.These criteria come with some strings attached, as they can contaminate the solution.A broad family of such adaptive criteria have been formulated using the finite-volume (FV) [9,10,14,15,17,[22][23][24][25][26][27][28][29][30] , the discontinuous Galerkin (DG) [31][32][33][34][35] , and the Flux Reconstruction (FR) [36,37] mathematical frameworks.One of the more promising techniques that can be used across any of the three numerical approaches is the MOOD technique, first conceived by Diot et al. [38] .The paradigm has since improved, see Diot, Clain, Loubére, Dumbser, Nogueira and others [22,[38][39][40][41][42][43][44][45] and found to be successful in a variety of difficult flow configurations.An important ingredient is that the candidate solution at any cell is subjected to several test criteria; only those solutions which satisfy these criteria are admissible.Conversely, the cells with non-admissible solutions have to revert to either a lower-order approximation or another type of scheme or either another numerical framework (hybrid DG-FV implementations) until these criteria/conditions are satisfied.The criteria/conditions take jointly into account the physical admissible solutions and the numerical admissible solutions, and the reader is pointed to previous work for a comprehensive analysis on the formulation of these criteria [38,45] .
The most prominent drawback of the MOOD paradigm in its original incarnation is the associated computational overhead once cells exist which do not satisfy the predefined criteria.The overhead emerges first from the requirement to revert to another method/scheme/order for the afflicted cells.It is not, however, contained in those cells, and becomes a very pronounced issue since the solutions predicted at the neighbours of this cell has been contaminated through the flux reconstruction.Therefore, the ultimate goal of the numericist is to tag as fewer cells as possible without degrading the robustness of the framework; a challenging compromise.For this reason, the starting numerical method is going to be the Compact Weighted Essentially Non-Oscillatory (CWENO) scheme, which belongs to a type of methods for unstructured meshes introduced by Tsoutsanis and Dumbser [9] , that have been extended to multicomponent flows [10] .The non-oscillatory properties of this type of method will result in a reduction of the number of cells that do not satisfy the MOOD criteria at the outset, as opposed to an unlimited high-order method.Other notable numerical methods that could have been used include the family of Target Essentially Non-Oscillatory (TENO) methods that have been successfully applied to a wide-range of prob-lems [46][47][48][49][50] .The augmented robustness provided by the MOOD technique is necessary for flow problems with very strong gradients that are common at interfaces between solids, liquids and gases through which shock waves propagate.Every schemes and paradigm used is this work has been implemented by the community of the open-source solver UCNS3D [51] .This work will be assessing the performance characteristics of the code in a series of demanding 2-D and 3-D numerical tests.
The paper is structured as follows.In Section 2 we will be describing the formulation for the governing equations we opted for this study, after which we introduce the numerics used, the selected discretisation options for the fluxes and temporal advancement, and the MOOD algorithm in Section 3 .In Section 4 the numerical results obainned from the selected computational tests will be discussed and evaluated using available analytical results, reference solutions or experimental data for our test problems.In closing, our conclusions will be summarised in Section 5 .

Governing equations
For this study we employ the five-equation model of Allaire [8] which consists of the mass conservation equation of two fluid components continuity Eqs. ( 1) and ( 2) , the momentum equation and energy equation for their mixture Eqs. ( 3) and ( 4) respectively, and the volume fraction advection for one fluid Eq. ( 5) which is non-conservative as given below: ∂ ρu ∂t With ρ, p, u = (u, v , w ) T , E and a being the density, pressure, velocity, total energy and the volume fraction respectively.For closing the five-equation model we employ the stiffened gas EOS where the pressure for each state is: with the reference pressure being denoted as π ∞ ,i ≥ 0 , and for gases π ∞ = 0 .The total mass is given by: and ρ being where (internal energy) being defined as ρ E − 1 2 ρuu .The EOS of the mixture is provided by the following expression We adopt the approach of Johnsen and Colonius [17] for reformulating the volume fraction advection Eq. ( 5) as follows: 3. Numerics

The unstructured framework of UCNS3D
Considering a 3-D domain that is made up of any combination of element shapes that can include hexahedral, tetrahedral, prismatic and pyramidal elements where each one is uniquely identified by index i , the equations are written in the following form: where W = W (x , t ) represent the conserved variable vector and volume fraction, and the non-linear flux normal to the cells faces being F n as given below: where u n is the velocity normal to the cell-interface and the source term s contains only the contribution from the volume fraction term a 1 ∇ • u of Eq. ( 5) .We employ a surface integral approximation of the source term similarly to Johnsen and Colonius [17] that make use of the approximate Riemann solver velocity estimate (u n ) Riem. as follows: Adopting a high-order finite-volume formulation for each mesh element and integrating Eq. ( 13) we obtain the following expression: where W i is the vector of the mean value of the variables and F n io is the numerical flux across the cell interface normal direction between (considered) cell i and neighbouring cell o and a n i, 1 being the mean volume fraction for cell i .N f corresponds to the number of element faces, N qp the number of surface quadrature abscissas, and | S io | the area of the considered face.The high-order solutions at the cell interfaces for cells i and j are given by W n io,L (x io,α , t ) and W n io,R (x io,α , t ) respectively.A Gauss-Legendre quadrature is used with α denoting the abscissas x α and weights ω α for each cell interface.

Spatial discretisation
In the current study with UCNS3D [51,52] the least-square reconstruction approaches of Tsoutsanis et al. [9,10,23,53] and Titarev et al. [54] are adopted, that have been applied to laminar, transitional, and fully turbulent flows at the a wide-range of Mach numbers [9,10,23,45, .Herein the most key details of this approach are presented, and the reader is pointed to previous work for a detailed description of the implementations For an arbitrary cell i obtaining a high-order r + 1 order, requires the formation of an r order polynomial p i (x, y, z) which has to provide the volume-averaged value W i as follows: Since we are dealing with unstructured meshes that can contain elements of different shapes and sizes, the transformation from the physical space (x, y, z) to a reference space (ξ , η, ζ ) can reduce the influence of the scaling effects introduced by different shapes and sizes of elements.Each of the elements that is not a triangle or tetrahedral is decomposed into them just for the reconstruction process [53] since arbitrary unstructured hexahedrals, prism and pyramids might not be transformed to their unit-length counterparts as detailed in Tsoutsanis et al. [53] .During this transformation the volume average of a variable remains constant as indicated below: wit i being the considered cell's volume in (ξ , η, ζ ) space.For obtaining the coefficients of the polynomial we make use of the information contained in the close spatial neighbourhood of the considered cell by forming a central stencil S 1 , where it is formed by adding M + 1 neighbouring elements.For an overview of different stencil selection strategies the reader is pointed to the work of Tsoutsanis [69] .Herein we employ a stencil-based compact algorithm (SBC) from Tsoutsanis [69] .
The number of the unknown coefficients of the polynomial are denoted by K and given by: with d ∈ [2 , 3] denoting the space dimensions.A union of all the M + 1 elements composes the central stencil S 1 as follows: where m is corresponding to the numbering (locally) of the stencil elements, 0 being reserved for the first cell which always the one considered, and c = 1 referring to the central stencil.For avoiding ill-conditioned linear systems for obtaining the polynomial coefficients we use twice as many elements in the stencil as polynomial coefficients M = 2 K, to ensure that an overdetermined system is formed which has been documented to improve the robustness of the reconstruction [23,24,69] .
The stencil is transferred to the reference space, and becomes S c i , and the polynomial can is provided by expanding locally the basis functions as follows: The vector of the conserved variables for cell i is W 0 and a k denoting the unknown polynomial coefficients (degrees of freedom).For each cell a k must adhere to the aforementioned condition in which the reconstruction polynomial (volume Therefore the basis functions ψ k are formulated in a way to automatically satisfy the Eq. ( 18) and they are defined as follows: In the current study the polynomials ψ k correspond to the Legendre polynomials.Defining the integral of the basis function as A mk over a cell m , and the right hand side vector by b as follows: (25) and rewriting them in a matrix form we get: The overdetermined linear-system is solved with the QR technique and pseudo inverse (Moore-Penrose) of A mk is prestored for computational efficiency as reported in Tsoutsanis [69] .

The MUSCL framework
The MUSCL class of schemes is used as one of the methods that the MOOD algorithm can resort to when other higherorder methods are not working.For this class of methods all the polynomial coefficients are limited based on the violation or not of certain bounds that the reconstructed solution will need to be within.The scheme can be formulated as: (27) with W l,α being the projected reconstructed solution at cell interface l and quadrature (ξ a , η a , ζ a ) , with the volume average conserved variable being W i .The degrees of freedom are limited by θ i which is a unique limiter cell i .The MUSCL algorithm works in the following fashion: 1. Obtain the minimum and maximum bounds of each variable from cell i and its neighbours 2. For every face and quadrature point obtain the high-order reconstructed solution W l,α for every face and quadrature point 3.For every face and quadrature point obtain the maximum allowed limiter θ l,α 4. Use the minimum value of the limiter θ i = min (θ l,a ) so that each cell has one unique limiter value.
The Barth and Jespersen [72] limiter is used where it can provide up to 2nd-order of spatial accuracy, while for higherorder of accuracy the Michalak and Ollivier-Gooch [73] type of limiters (MOG) and variants of them [65] can be used.The reader is pointed to Tsoutsanis [65] for description of high-order limiters on unstructured meshes.

The CWENO paradigm
The CWENO implementation formulated by Tsoutsanis and Dumbser [9] is adopted for this work.The CWENO scheme combines a high-order polynomial associated with the central stencil and with directional stencils of lower-order of accuracy and polynomials to build an optimal polynomial p opt .In this fashion when the solution is discontinuous the lower-order polynomials will have a greater influence in the reconstructed solution, while in regions of smooth solutions the high-order polynomial will have a greater influence and as a results the optimal polynomial will be recovered.The CWENO variant is preferred from the traditional WENO variants [23] in terms of computational cost (lower-order polynomials are cheaper to compute) and non-oscillatory behaviour (lower-order polynomials less oscillatory).The optimal polynomial is given by: (28) with s being the index of the stencil, with the central s = 1 , and the directional ones s = (2 , 3 , . . ., s t ) , and the linear weight defined by d s is that their sum is required to be 1.For computing the p 1 polynomial we subtract the low order polynomials from the optimal one: The resulting CWENO polynomial is obtained by a combining all the polynomials non-linearly : (30) with ω s being the non-linear weight for every polynomial, and when the solution is smooth ω s ≈ d s , hence obtaining the high-order accuracy, when the solution is discontinuous a low-order solution dominated by the directional stencils will be obtained.The non-linear weight ω s are defined as: For this study we use = 10 −14 and b = 4 .The directional stencils use a linear polynomial that results in 2nd order while any desired order can be obtained from the central stencil.For a comprehensive documentation of the smoothness indicators for the CWENO scheme the reader is pointed to Tsoutsanis and Dumbser [9] .For the high-order central stencil λ 1 is given a large value and then all the weights are normalised (to ensure that their sum is unity) as follows: where all of the linear weights of the directional stencils are equal: where the number of the stencils is denoted by s t .

Fluxes & time advancement
For the numerical fluxes we employ the HLLC approximate Riemann solver introduced by Toro [74] , since it is well adapted for the volume fraction treatment employed in this work, and the for time advancement we use the 3rd-order Runge-Kutta method [75] , where a CF L = 0 .6 is used for all tests unless otherwise specified.The reconstruction is performed with respect to the primitive variables since they were previously found [10] to be less oscillatory than the conserved variables, and finally a suitable Gaussian quadrature rule for the chosen polynomial order is used for integrals approximation.

MOOD
The present study uses the multi-order optimum detection (MOOD) algorithm [45] .The algorithm evaluates a posteriori the solution of the high-order numerical method using a set of criteria that detect a variety of even minute oscillations.It is therefore possible to guaranty the solution remains monotonic and its overall physicality is also ensured.In all, the algorithm enhances the non-oscillatory properties of any numerical method.Previous work by Farmakis et al [45] resulted in a revision of this a posteriori paradigm where the admissible criteria are strategically relaxed to increase the performance of the code without altering the ability of the algorithm to capture and control spurious oscillations.Our relaxed-MOOD implementation has been extensively tested in demanding computational problems like the non-linear growth of the Rayleigh-Taylor and Richtmyer-Meshkov instabilities with sharp discontinuities and steep density, temperature and velocity gradients.The main platforms have been high-and intermediate-Atwood number supersonic 2-D and 3-D converging flows relevant to highenergy-density systems.Since these problems where done using a single γ value for the different parts of the flow, it has been a major objective to extend the framework to multicomponent flow problems.
The algorithm is special in the sense that all cells in the computing volume will need to have a candidate solution which is to be concurrently physically and numerically "admissible" before the solution may be updated.If either category of these conditions is violated in any cell, then the triggered sensor(s) call(s) up the MOOD algorithm to temper the solution according to a prescribed strategy, first outlined by Diot et al. in Clain et al. [38] , and then in Farmakis et al. [45] .The order of the solution's accuracy is reduced locally using the so-called Discrete Maximum Principle (DMP) -originally; it was subsequently optimized in Farmakis et al. [45] -until the candidate solution is admissible.This involves switching efficiently the numerical schemes and their order appropriately so that a minimal amount of numerical diffusion is injected into the target region, thus maintaining the global intended high-order of accuracy in the solution.To be more specific, every highorder accurate solution predicted in the computational domain at a given timestep by schemes like CWENO, MUSCL, etc (i.e. the candidate solutions) will have to be scrutinized by two detectors right after the solution of the fluxes and before solution update: 1. Physical Admissible Detector (PAD): The first detector, which is responsible for ensuring the physicality of the candidate solution.The detector will remove every negative pressure and density value in the computational domain, since otherwise the solution predict other non-physical sound speeds, imaginary time steps and so on.The physicality here is assessed from the point of view of a fluid flow, which limits the generality of the criteria as far as predicted pressures are concerned.2. Numerical Admissible Detector (NAD): Our second sensor correspond to a more versatile variant of the DMP (see Clain et al. [38] for more on the topic).This sensor is responsible for guaranteeing the monotonicity of the solution; that is, no spurious minima or maxima are introduced locally by the reconstructed solution.It compares this with the reconstruction estimated in the Runge-Kutta step prior, a crucial element of the implementation which ensures optimal operation of the method [45] .
If spurious oscillations have contaminated a candidate solution, it will fail to pass the PAD &/or the NAD criteria.As a consequence, the MOOD algorithm drives the code to locally downgrade the order of accuracy by using an auxiliary scheme.The implemented method uses the 2nd-order MUSCL scheme for its auxiliary numerical method.Given that the oscillations may persist, the PAD & NAD criteria are applied again to the new downgraded candidate solution.In a crucial departure over the original proposal by Diot et al. [38] , UCNS3D 's version [45] apply now both kinds of admission criteria using relaxed tolerances.Our solver's variation, using more relaxed criteria has been carefully chosen so that the code can simultaneously ensure that the predicted solution in a given cell is a admissible and has the highest order of accuracy.If the solution satisfies the new criteria, then it stabilized and may move on in time.Otherwise the solution failed to satisfy the relaxed criteria and the code parachutes to the bulletproof scheme.For UCNS3D , this scheme has been selected to be the Godunov's upwind scheme; it is first-order accurate and satisfies both admissible criteria criteria by default; the solution can advance in time.While, in theory, the algorithm can propagate the solution though a chain of multiple different schemes and orders of accuracy, it was found during implementation by Farmakis et al. that the cascade should be maintained as short as possible, otherwise the footprint of the method offsets any benefits in stability.We should point out that MOOD's a posteriori nature necessitates the reevaluation of the solutions for the failed cell's neighbours.This is so because their solutions can still be influence by the troubled cells due to the fluxes, and carry this information in future time steps, hampering the effectiveness of the method.
The algorithm of the relaxed-MOOD paradigm implemented by Farmakis et al. may be found in Fig. 1 .In this flow chart, the group of cells V i represents the von Neumann neighbours for the stricken cell i targeted by the algorithm.In the current formulation, the NAD criterion evaluates the vectors of the conservative variables, following the example of Diot et al. [22] .During every Runge-Kutta stage temporally advancing the solution, it can be seen that the candidate solution U * i at the targeted element i must be bounded by a certain range of values.This is provided by the V i region heretofore defined.It should be clarified at this point that the superscript ( n ) in the chart defines the previous Runge-Kutta stage, rather than the previous time step: min Fig. 1.MOOD-algorithm implementation flow-chart for UCNS3D code.The user-selected high-order numerical method (such as CWENOZ, WENO etc.) passes on the information of the fluxes to the a posteriori MOOD algorithm.The switch between criteria for the PAD/NAD sensors is shown in steps 5, 6a for higher than 2nd-order methods, and 6b for 2nd-order method.Transitioning from a high order method to the first auxiliary (2nd order MUSCL) occurs in step 9a, while at step 9b the solution resorts to the first order Godunov flux.If the first order method the PAD/NAD checking is no longer deployed.
Table 1 e L ∞ and e L 2 errors, and convergence orders of the schemes with respect to volume fraction and normalised wall-clock time for the 2-D multi-species convergence test after one period.It can be noticed that the MOOD-CWENO variant achieves the same order of accuracy as the CWENO variant, with a small computational overhead (3-12%) due to the PAD/NAD condition checking.The relaxed DMP margins, used by the implemented scheme to reduce the order of accuracy from the originally used to second-order can be seen in step 6a. of Fig. 1 , marked as ( • ) :

CWENO MOOD-CWENO
whereas the more lenient criteria, employed in step 6b. of Fig. 1 and marked by R , allow a smooth transition along the scheme/order of accuracy cascade without excessive introduction of numerical dissipation by marking as fewer cells as possible.In the current implementation these lenient criteria are used when the solution is downgraded from second to first order of accuracy: Using the relaxed criteria in UCNS3D made it possible to drastically reduce the computational footprint of the MOOD method, while at the same time leaving the framework's effectiveness unaffected.In conjunction with a desireable nonoscillatory behaviour of the used limiters by the MUSCL scheme [65] , it was shown during testing of UCNS3D 's MOOD implementation in Farmakis et al. [45] that the algorithm guarantees the non-oscillatory properties anticipated by the solution.
It has to be stressed that for the present MOOD implementation only the density, total energy and volume fraction are subject to the NAD criteria.For a comprehensive description of the implementation the readers are pointed to the work of Farmakis et al. [45] .One has to be aware that the stiffened gas (EOS) chosen might admit solutions with negative pressure.The MOOD algorithm might tag cells because of its intrinsic PAD conditions even though the solution is valid under the stiffened gas EOS.While, generally, predicting negative pressures in the solution implies that the material is a solid, this is problematic for a fluids code since it would result in imaginary sound speeds blowing up the solution [76] .It is important that the user is cognizant of the nuances that come with the EOS employed each time, as it will be seen in some of the results presented in this work.

Multicomponent convergence test
The primary goal of this test is to understand the behaviour and influence of the MOOD technique for a flow problem without sharp gradients and demonstrate whether it can achieve the designed order of spatial accuracy of the selected numerical method.We are also interested to understand the computational cost associated with the conditions/sensors/criteria checking of the MOOD variant.The multicomponent advection of a smooth volume fraction of two gases (nitrogen and helium with specific heats γ = 1 .4 and γ = 1 .66 ) for one period is considered as introduced by Wong and Lele [16] .The initial condition is given by: The 2-D computational domain [0 , 1] 2 is discretised by an unstructured triangular mesh with a resolution of an approximate edges length h of 0 . 1 , 0 .05 , 0 .025 , and0 .0125 , and the initial profile is advected for one period which corresponds to a time of t f = 1 .0 (normalised).The numerical errors in terms of the e L 2 and the e L ∞ are examined with the computed solution being    W c x, t f and the exact solution being W e x, t f .The initial condition is the exact solution W e x, t f since the initial profile returns to its initial location after one period.Since we want to minimise the time-discretisation error and focus on the spatial discretisation, a CF L = 0 . 1 is used.For appreciating the computational cost the wall-clock time that each simulation took was normalised with respect to the cheapest numerical scheme (CWENO3) on the coarsest mesh.
As expected the MOOD-variants achieve the same errors and order of convergence as the CWENO schemes for this smooth test problem as seen in Table 1 .This is because our MOOD algorithm remains dormant, since no cells have nonadmissible solutions.However just the checking the admissibility of the solutions translated to a small computational overhead (3-12%) due to the PAD/NAD condition checking.

Helium bubble -shock wave interaction
The widely-used test problem of a helium bubble interaction with a shockwave which was introduced experimentally by Haas and Sturtevant [77] is considered herein in 2D.It is a popular test problem [14,15,17] for validation and verification of numerical methods, models and techniques for multicomponent compressible flows.The arrangement involves a shock tube filled with air ( γ = 1 .4 ), containing a helium ( γ = 1 .66 ) bubble (not pure but with 28% mass concentration of air) of diameter D b = 5 cm that is hit with a shockwave as shown in Fig. 2 .The initial conditions are given by: for Pre-shock ( 0 . 0 , 1. 658 ,−114 . 49 , 0 , 159060 , 0) , for Post-shock ( 0 . 158 , 0. 061 , 0 , 0 , 101325 , 0. 95) , for Bubble.(38) A hybrid unstructured mesh is used that has triangles and quadrilateral elements with an average cell edge length of e m ≈ D b / 120 at the centre of the domain where the shock-bubble interaction takes place.For the right and left boundary a supersonic inflow and outflow is enforced, and for the rest of the domain reflecting boundary conditions are applied and the problem is simulated for t = 10 0 0 μs.The primary goal of this particular physical configuration is to evaluate the influence of the parameter used for the second cascade of the NAD sensor and in particular we explored three values (0.1, 0.05 and 0.5) for what is going to be referred as MOOD1, MOOD2, MOOD3, respectively.This parameter's value is applicable only to the NAD sensor for the second cascade from 2nd-order MUSCL to first-order.Thus, the objective here is to demonstrate how sensitive the results are with respect to this parameter, since it has a direct bearing on the number of cells that might end-up being solved by a first-order method.The selected high-order scheme is the 5-order accurate CWENO method.
The simulation results from all the methods exhibit the expected flow structures associated with this problem such as the vortex ring, the Kelvin-Helmholtz instabilities and the jet as shown in Fig. 3 which agree well with other studies [10,14,15,17] .What is important here is that the MOOD2 (which has the most-sensitive NAD sensor) activates more cells for first-order as seen in the contour plots of the instantaneous order of accuracy, whereas the other methods have a negligible number of cells solved by a first-order method.From the position-time histories of the interfaces as defined in Fig. 4 it can be noticed that there are not significant differences between the MOOD schemes, and are in close agreement with several reference results [10,78,79] .Fig. 6.Contour plots of density gradient magnitude (left), volume fraction(middle) and mood enabled regions (right) of the shock wave interaction with a water column test, at various instants (left to right).It can be seen that some regions have resorted to the first-order bulletproof scheme primary due to the presence of negative pressures admitted by the selected EOS.[80] , where a good agreement is achieved for the single instance ( t * = 0 .8 ) results provided.

Water column in air
This test involves a shock wave interacting with a cylindrical water column that has been investigated computationally [10,80] and experimentally [81] .We are interested to understand how the MOOD algorithm performs in this test that involve very strong gradients sicnce a shockwave M sh = 2 .4 in air hits a water bubble.The initial state of the test problem can be seen in Fig. 5 .
The initial condition is given by: 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 10 0 0 , 0 , 0 , 10 5 , 6 .12 , 0 .343 • 10 9 , 1 , for Water (39) The domain is x ∈ [0 , 0 .2208] , y ∈ [0 , 0 .1152] m is discretised with a quadrilateral mesh of 891,497 cells, with nonreflecting boundaries at the sides and reflecting boundaries for the rest of the domain.The water column is placed at (0 .0576 m , 0 .0576 m ) , and the shocked region is between 0 and 0.05 m with the rest of the domain being at Pre-shock conditions.For this problem the CWENO3-MOOD scheme is used and and the flow problem is simulated for t * = tu / D b = 10 , with u being the velocity initially at the pre-shocked region.From the results obtained as shown in Fig. 6 the typical flow structures expected are well resolved and they include the transmitted wave, the expansion wave and the Mach stem [10,80,81] .
It can be seen in Fig. 6 that some cells are tagged for solution by a first-order method, but this is not due to the NAD sensor, but primarily due to the fact that the stiffened gas EOS can allow solutions with negative pressure.The EOS is responsible for inadvertently triggering the PAD sensor (since negative density or pressure are not allowed due to the fundamental assumption of a fluid flow in the current MOOD paradigm) and as a result, the solution resorts to the first-order method.It is obvious that the appearance of the negative pressure is not resolved by the bulletproof first-order methods, since it is a physical prediction of the EOS itself .
Regardless, it is also evident by the comparison of the present work with [80] (see Fig. 7 ) that we obtain valid results.The computed density gradient magnitudes are in good agreement with the numerical prediction in Xiang and Wang [80] (this may be found in Fig. 6 ).We note that the figure shows that the MOOD algorithm switched some regions to the first-order bulletproof scheme.As we explained previously, this is primarily due to the presence of negative pressures admitted by the selected EOS.As also documented by previous computational and experimental studies [10,[80][81][82] our results exhibit the flow separation, the compression and flattening of the water column at late times.

2D underwater explosion
A two dimensional underwater explosion is considered as described in Nourgaliev et al. [83] , Hu et al. [84] , Fu et al. [85] .This problem involves an explosive gas cylindrical area of radius r = 0 .12 m placed underwater, that moves the water-air interface.The initial condition is given by: ( ρ, u, v , p, γ , π ∞ , a 1 ) = ( 1 . 225 , 0 , 0 , 101325 , 1. 4 , 0 , 1) , for air 1250 , 0 , 0 , 10 9 , 1 .4 , 0 , 1 , for shocked gas ( 10 0 0 , 0 , 0 , 101325 , 4 . 4 , 0) , for Water (40) We are interested to understand how the MOOD algorithm performs in this challenging test problem.For this reason we have employed a CWENO5 method, and the MOOD algorithm is modified to neglect regions with negative pressure that originate from the stiffened gas EOS, therefore making the algorithm more relaxed.Looking at the obtained results in Fig. 8 it can be noticed firstly that the correct flow pattern is obtained which is close agreement with the results of Nourgaliev et al. [83] , Hu et al. [84] , Fu et al. [85] .It has to be stressed that the cells that are tagged by MOOD are not significant, and more importantly the majority of the tagged ones are solved by a 2nd-order MUSCL scheme which is sufficient for satisfying the MOOD criteria.

Shock-induced collapse of 2-D bubble arrays
We extended our numerical studies by investigating the collapse of an array of bubbles initiated by an advancing shockwave.The goal of this test case is to understand the performance of the developed MOOD-CWENO algorithm in this more complicated physical configuration.The main element of interest is the strong accumulation of energy as the shockwave implodes this constellation of bubbles.We are considering the 2-D arrangement of an array of three bubbles arranged in a triangular shape, as introduced by Bempedelis and Ventikos [86] , shown in ( ρ, u, v , p, γ , π ∞ , a 1 ) = ( 1 , 0 , 0 , 10 0 0 0 0 , 1 . 4 , 0 , 0) , for gas bubbles 1323 , 0 , 0 , 10 0 0 0 0 , 4 .4 , 6 • 10 8 , 1 , for Water (41) The domain is discretised by three quadrilateral meshes of resolution equivalent to 20 0, 10 0, and 50 points per large bubble radius respectively.The simulation is run for t = 1 .2 μs with a CWENO5-MOOD method.From our simulation results, at different times and for the finest mesh, we show in Fig. 10 that our results demonstrate similar trends to the results of Bempedelis and Ventikos [86] (even if the Mach number is slightly different for the present study) in terms of the predicted flow-structures.The maximum pressure histories from our simulations are shown in Fig. 11 , where different grid resolutions are used.For the finest mesh, a maximum pressure peak of approximately ×6 .3 the pressure at the shocked region is achieved close to the centre of the small bubble at approximately t = 0 .93 μs, shown in Fig. 10 .In the same figure, it can be noticed that a significant portion of the domain is solved with a first-order method, and, again, this is partially attributed to the negative pressure allowed by the selected EOS.Therefore the PAD sensor inappropriately tags these cells to successively reduce to MUSCL and the first-order method.This "overcorrection" cannot resolve the otherwisephysical prediction of the stiffened-gas EOS.This implies that it might be worth redefining the PAD sensor to remove any such overreaction for a selected EOS.

Shock-induced collapse of 3-D bubble arrays
The final numerical test case is an extension of the previous numerical study, to the third dimension, and a parametric investigation of different shapes, sizes and arrangements of the gas regions.The bubble configurations considered in this series of tests are the Triangular (Triangular-Loose), the Pyramidal, and the Toroidal (Toroid-Loose) arrangement, as investigated by Bempedelis and Ventikos [86] .In addition to these, we expanded our parameter space and investigated certain variants of the original test cases; we developed a variant of the Triangular arrangement, which we term as Triangular-Compact, a modified Toroidal configuration (Toroid-Compact) and, finally, we experimented with an Ellipsoidal bubble.The   initial profiles for all different bubble shapes and configurations are presented in Fig. 14 , and the corresponding parameters are given in Table 2 .
In the previous Section 4.5 , the regions of the air cavities/bubbles have been generated using the equation for area within a circle,

Table 2
Parameters of the 3-D bubble configurations.obtained by Bempedelis and Ventikos [86] , as presented in Fig. 15 (neglecting the mismatch in Mach numbers between the present study and [86] ).
The results using the finest prismatic mesh in the ellipsoidal configuration are shown in Fig. 16 .The number of cells that are reconstructed using a first-order method is negligible compared to the previous 2-D configurations, and most of the cells are solved by the CWENO4 variant employed for this 3-D problem.For this problem we have chosen to decouple the negative pressure values from the PAD sensor, since the selected EOS is the one responsible for them.As a result, the simulations were performed without any PAD-wise over-correction, and the NAD sensor was the criterion mostly dictating the behaviour of the MOOD algorithm.

Conclusions
This paper explores the extended applicability of the MOOD algorithm to high-order CWENO schemes on unstructured meshes for compressible multicomponent flows.The desirable features of the CWENO schemes associated with robustness and accuracy are well suited for multicomponent flows since fine flow structures are well resolved and the solutions are essentially free from spurious oscillations.The fortification of the selected numerical method with the MOOD implementation benefits the robustness of the numerical framework making simulations of stringent flow problems feasible.The MOODenabled CWENO schemes can deal with the very strong gradients encountered in some of the test problems investigated such as the of bubble array collapse.Especially in the 3-D shock-induced collapse of the bubble arrays the most stiff of the problems explored, the framework showed unwavering robustness.It has also been found that in the presence of an EOS that physically admits negative pressures, the PAD sensors attempt to over-correct the solution.This points out that the existing PAD formulation should be redefined in physically motivated ways which will ensure that the EOS predictions are correctly accounted for in the solution while the entire MOOD implementation maintains its demonstrated robustness.Future development explore the application of the MOOD-CWENO schemes to more complicated flow problems in 3-D using further enhanced diffused interface models, and redefined admissible solution sensors to verify the generality of the MOOD implementation.

Fig. 2 .
Fig. 2. Computational domain setup for the helium bubble interaction with a shock wave test.

Fig. 3 .
Fig. 3. Time histories of the interfaces position obtained with CWENO5-MOOD schemes and comparison with the computational results from Terashima and Tryggvason[78] , Quirk and Karni[79] , and Tsoutsanis et al.[10] .All MOOD variants agree well with the reference results without significant differences between them.

Fig. 4 .
Fig.4.Plots of the volume fraction contour (top) for the interaction of the helium bubble with a shock wave test at t = 983 μs , and the MOOD tagged cells and their corresponding attainable spatial order for the MOOD1 (left), MOOD2 (middle) and MOOD3 (right) variants.It is clear the MOOD2 variant has an increased number of cells with 1st-order of accuracy as expected, while the other two have negligible number of cells, while not significant differences can be noticed between them regarding their solution profiles.

Fig. 5 .
Fig. 5. Drawing of the initial configuration for the shock wave interaction with the cylindrical water column test.

Fig. 7 .
Fig. 7. Plot of the pressure variation along the centreline of the water column in air flow problem, at different times, and comparison with the results obtained from Xiang and Wang[80] , where a good agreement is achieved for the single instance ( t * = 0 .8 ) results provided.

Fig. 9 .
Fig. 9. Schematic diagram showing the setup for the triangular 2-D bubble array test problem.

Fig. 9 .
The domain extent for the test problem is x ∈ [0 , 0 .003] m , y ∈ [0 , 0 .003] m .The radius of the small bubble is R b = 0 .0002 m and the radius of each greater bubble is set to R B = 0 .0 0 05 m .Post-shock conditions correspond to a Mach number of M = 1 .72 , and the initial conditions for both gases read:

Fig. 10 .
Fig. 10.Contour plots of density (left), pressure gradient (middle) and MOOD activated cells (right) for the 2-D bubble array shock-induced collapse test problem at different times for the finest mesh.It can be seen that the maximum pressure is achieved close to the centre of the small bubble at approximately t = 0 .93 μ, and that several cells have been solved by a first-order method.

Fig. 11 .
Fig. 11.Time history of maximum pressure in the domain for the bubble array shock-induced collapse test problem at different grid resolutions.

Fig. 12 .
Fig. 12. Cutaways of the fine prismatic mesh used for the shock-induced collapse of 3-D bubble arrays.

Fig. 13 .
Fig. 13.Evolution of the maximum recorded pressure in the computational domain for the 3-D bubble array shock-induced collapse test; Different configurations on the coarse hexahedral mesh (left), and Ellipsoidal configuration on the prismatic fine and coarse hexahedral mesh (right).

Fig. 14 .
Fig. 14.Schematic diagrams of the 3-D bubble configurations investigated, from top (left to right); Triangular-Loose, Triangular-Compact, Toroid-Loose, Toroid-Compact, Ellipsoidal, and Pyramidal; v1 is the density gradient magnitude, v2 is the pressure gradient magnitude and the iso surfaces of the volume fraction are coloured by the volume fraction.

Fig. 15 .
Fig. 15.Contour plots of pressure ( yx -plane), volume Fraction ( yz -plane) and iso surfaces of volume fraction flooded with density at different times for the pyramidal configuration.

Fig. 16 .
Fig. 16.Contour plots of pressure (left-plane) at the middle of the domain, density (right-plane) at the middle of the domain, MOOD activated cells (bottom plane) at y = 0 .002 , and volume fraction coloured by pressure for the ellipsoidal configuration of the shock-induced collapse test problem at different times.