Cavitation Bubble Cloud Break-Off Mechanisms at Micro-Channels

This paper provides new physical insight into the coupling between flow dynamics and cavitation bubble cloud behaviour at conditions relevant to both cavitation inception and the more complex phenomenon of flow “choking” using a multiphase compressible framework. Understanding the cavitation bubble cloud process and the parameters that determine its break-off frequency is important for control of phenomena such as structure vibration and erosion. Initially, the role of the pressure waves in the flow development is investigated. We highlight the differences between “physical” and “artificial” numerical waves by comparing cases with different boundary and differencing schemes. We analyse in detail the prediction of the coupling of flow and cavitation dynamics in a micro-channel 20 μm high containing Diesel at pressure differences 7 MPa and 8.5 MPa, corresponding to cavitation inception and "choking" conditions respectively. The results have a very good agreement with experimental data and demonstrate that pressure wave dynamics, rather than the “re-entrant jet dynamics” suggested by previous studies, determine the characteristics of the bubble cloud dynamics under “choking” conditions.


Introduction
Cavitation has been the subject of intensive research for decades due to its relevance to many engineering applications such as micro-channel flows, fuel injectors and propellants [1,2]. Initial stages of cavitation consist of the generation of a cavitation cloud attached to the wall that primarily develops within the shear layer of the flow and is largely unobtrusive to the main flow [1]. As the cavitation region evolves, and once a critical cavitation number is reached, cloud break-off occurs. During this process the cavity volume oscillates in either an intermittent or periodic manner. Understanding the mechanisms that govern the cloud break-off process, as well as the frequencies at which it occurs, is important since cloud break-off is linked to effects such as loss of lift in hydrofoil applications, and more generally to unwanted vibration and erosion. We recall here that the cavitation number (Cn) characterises the onset and extent of cavitation and is the ratio between a reference pressure (P 0 ) and the vapour pressure (P v ), P 0 − P v , divided by the dynamic pressure of the free stream flow, ρv 2 ∞ /2. Lower values indicate a higher probability of cavitation [3]: The occurrence of cloud cavitation break-off has been related in many previous studies to the presence of a "re-entrant jet flow" at the downstream region of the cavity which forms as the liquid flow outside the cavity reattaches to the closure region [1]. One of the earliest studies that demonstrated this phenomenon is by Furness and Hutton (1975) [4]. They used a potential-flow-based numerical method in order to compute the structure of the re-entrant flow on a convergent-divergent nozzle until it intersected the cavity interface. Franc [1] described the existence of a re-entrant jet as a necessary outcome of the cavity region being the region of minimum pressure. Thus, the pressure gradient around the cavity region is directed away from its lowest point, forcing the flow along the liquid-vapour interface to curve towards the cavity. The jet originates from the cavity closure region and exists as a thin viscous film travelling upstream over a surface [5][6][7].
More recently, it was suggested that with a sufficient reduction in the cavitation number (giving rise to greater cavitation propensity), large-scale, periodic cloud shedding can be observed and associated with the formation and propagation of a shock within the high void-fraction bubbly mixture in the separated cavity flow [5,6,8,9]. Jahangir et al. [6] established within a venturi nozzle that cloud break-off was produced by the presence of the re-entrant jet for a cavitation number Cn > 0.95 while, for Cn < 0.75, cloud shedding was driven by bubbly shock effects generated by the collapse of previously shed cavities. This transitional regime was also demonstrated by Pelz et al. [7] where a relationship was established between the ratio of the Reynolds (Re) and the cavitation numbers, and the occurrence of this transitional regime. In the work of Ganesh et al. [5,9] with the use of high-speed imaging and time-resolved X-ray densitometry of cavitation over the apex of a wedge, the features of the observed condensation shocks were related to the average and dynamic pressure, and the void fraction. The transition to strongly shedding conditions was associated with a drop in the speed of sound in the bubbly mixture within the partial cavity. The strongest shock waves occurred for Mach numbers higher than one. More specifically, as shock fronts from collapsed cavities impinged on a growing vapour cloud, the internal pressure wave propagated at a speed much closer to, or higher than, the local speed of sound which caused an internal bubbly shock and drove cloud break-off.
Both mechanisms of cloud break-off (re-entrant jet and/or shock waves) are challenging to capture numerically. The difficulty in the numerical study of the re-entrant jet is that it is mostly a boundary layer phenomenon that requires high grid resolution close to the wall. The mechanism of cloud break-off based on shock wave propagation is also hard to predict since time steps of the order of nanoseconds are required, along with fully compressible codes, to capture local pressure peaks [10]. An additional difficulty is that there are only a few qualitative experimental data sets available for cavitating flows (especially for micro-channel flows) and even fewer ones that clearly identify the mechanisms for cloud break-off. An important reason for this is the difficulty in obtaining adequate visualisation. For example, the work of Ganesh et al. [5,8,9] mentioned above, was performed for cavitation over an apex of a wedge and thus it is indicative of the dynamics of the mechanism of the cloud break-off in external geometries, while the work of Jahangir et al. [6] was performed in a venturi geometry that is much larger than fuel injectors. An experimental investigation that was performed for numerical validation of Diesel fuel injection processes was conducted by Winklhofer et al. [11]. The study presented 'line of sight' distributions of the vapour void fraction, produced by shadowgraph imaging, in a plane micro-channel between an upstream and a downstream plenum chamber, with the flow 'sandwiched' between two windows. Measurements exist for several different pressure differences between the chambers and this case is commonly used for numerical validation due to the simple orthogonal geometry and realistic Diesel injection pressure conditions [12][13][14][15]. The limitations of the shadowgraph images though should be noted; the visual data provided are averaged over 20 different time-steps and are the 'line of sight' average over the depth of the cross-section. Mauger et al. [16] report similar visualisations of an orthogonal geometry, providing good instantaneous images of the cavitation distribution along the centre of the nozzle. The experiments show that instabilities in the shear layers intensified as the pressure difference increased and produced vortices in the flow. Mauger et al. also showed the presence of pressure waves due to the collapse of vapour bubbles in the throat as well as at the nozzle exit. The averaged cavitation topology seen in this investigation is similar to that of Winklhofer et al. [11]. It should be pointed out, though, that in neither of these two experimental studies the cloud break-off dynamics were investigated.
Apart from elucidating the mechanism of cloud break-off, another important aspect is the identification of the frequencies of the pressure and vapour fluctuations associated with this phenomenon. A number of papers have studied the link between frequencies of pressure fluctuations and cavitation dynamics. Experimentally, Crua and Heikal [17] measured the injection process in a Diesel injector, using 3D laser Doppler vibrometry, with a focus on the distribution of the flow dynamics over time. Some peak frequencies were identified to be dependent on the internal nozzle cavitation along with other peak frequencies that were dependent on the injector nozzle geometry. A more recent effort to visualise vortical cloud shedding structures within an internal nozzle set up was attempted by Mitroglou et al. [18] who showed that an increase in velocity, and hence an increase in Reynolds number (Re), led to lower cloud shedding frequencies with a Strouhal number (St) number of 3.1 at the dominant frequency. In the study of Mihatsch et al. [19], a numerical simulation using a density-based finite volume method, taking into account the compressibility of both phases and resolving collapse-induced pressure waves (CATUM code), was performed on an airfoil geometry described in [20]. In their study, four characteristic frequencies were identified at all operating conditions (instead of two in the case of Peters et al. [2] for the same geometry). The first and third frequency in order of amplitude were found to increase with increasing pressure level (lowest frequencies). The frequency peak with the second highest amplitude was the shedding frequency and the fourth highest peak was found to be independent of the operating condition and related to the reflection of collapse-induced pressure waves at the radial boundary of the gap (8-9 kHz). Experimental visualisation and LES modelling of a vertical orthogonal nozzle geometry was performed by He et al. [21]. The authors used FLUENT (pressure based code) as their multiphase solver and the Schnerr-Sauer cavitation model. Two frequency peaks were identified to be caused by cavitation. The lower peak was reported to be related to cloud shedding at 3 kHz and the higher peak was related to the collapse of small bubbles at 6 kHz.
The unsteady nature of cavitating flows (including the phenomena relevant to the cavitation cloud break-off at the macro-scale as well as bubble collapse at the micro-scale) generates pressure waves that travel throughout the domain [22,23]. These propagating waves will interact with the domain boundaries in different ways (varying from 100% reflection to 100% transmission) depending on the experimental set up. When performing a simulation, these waves also interact with the numerical domain outlet (which might have to be chosen deliberately to not be the same as the real geometry outlet) and may be artificially reflected back (artificial in the sense that they would be transmitted during an experiment). This is particularly the case when the full domain is not modelled to reduce computational cost. To the best of the authors' knowledge no previous studies have been performed on the role of the boundary conditions in the prediction of cavitation dynamics, especially in terms of cloud break-off, apart from the preliminary work of Pearce et al. [14]. In this work we adopted a non-reflective boundary condition developed by Poinsot and Lele [24] within a compressible multiphase solver. This boundary condition attempts to resolve a far field pressure boundary condition producing a pressure gradient based on the distance to the far-field. Without the non-reflective boundary condition, propagated pressure or shock waves produced by cavitation collapse would be reflected back toward nozzle exit. In order to highlight the sensitivity of cavitation cloud-off to reflected pressure waves, a fixed value pressure condition was also adopted in our work alongside the non-reflective condition.
Apart from the boundary conditions, another consideration with the numerical representation of the propagating waves in cavitating flows is the damping effect introduced by various differencing schemes. The first order upwind scheme offers a fully bounded solution but is far too diffusive, and the second order central scheme has better accuracy but is unbounded. The central differencing scheme was used by Magagnato and Dumond [25] to simulate cavitation within a number of different geometries, producing a flat cloud topology with re-entrant driven cloud break-off. The upwind scheme, which is first-order accurate, was selected by Salvador et al. [26] to model cavitation in the experiment by Winklhofer et al. [11] using the Homogeneous Equilibrium Model (HEM). This scheme was selected for its stability but also produced a cloud topology which was flat and parallel to the wall with no convergence towards the centre of the duct as is seen in the experiment. Chen and Oevermann [27] modelled cavitation in the same geometry using a compressible multiphase stochastic framework. The differencing scheme used for the convective terms was a second-order upwind scheme which also produced a flat cavitation cloud topology and did not show any waves produced by cavity collapse. He et al. [21] conducted simulations of the same geometry using a compressible code with the Homogeneous Relaxation Model (HRM) to model phase change. A first-order upwind scheme was used for this investigation with the cavitation cloud appearing only within the laminar layer. Kärrholm, because of this, used a Total Variable Diminishing (TVD) scheme called MUSCL in OpenFOAM (taken from the name of the general sub-class of TVD solvers, monotone upstream-centered schemes for conservation laws) which maintains the monotonicity of the flow and thus is better at reproducing local large gradients [12]. These schemes resolve a blended scheme with a dynamic face limiter to suppress the higher order contribution at points of instability. Monotonic TVD schemes can handle discontinuities that arise from discontinuous functions or numerical inconsistencies [28]. Kärrholm's results show a cavitation cloud topology that was qualitatively more representative of the experimental visualisations and produced pressure and velocity profiles that were in good agreement with experiment. A three phase compressible solver with phase change across a single interface was formulated by Yu et al. [15] who conducted internal nozzle cavitation simulations using the Winklhofer geometry as well. A Normalised Variable (NV) gamma scheme was selected for the discretisation of the volume fraction and produced centreline pressure profile in very good agreement with experiments. The cloud topology though, appeared somewhat artificial. The analysis was not extended to cloud shedding. A parametric study of a high-resolution NV differencing scheme was performed by Guedri et al. [29] and found that, although the gamma scheme was computationally efficient, it was less accurate than other higher-order schemes such as MUSCL. Higher order schemes such as the Weight Essentially Non-Oscillatory (WENO) scheme are able to capture the large discontinuities and smooth regions [30]. However, it has been shown that higher order schemes require high computation and memory demand that can lead to stability problems [31].
Based on the previous observations, in the first half of our paper we present a sensitivity analysis in order to further examine the effects that the boundary conditions and the differencing schemes have on the predictions of the cavitation dynamics in a micro-channel with emphasis on the cavitation bubble cloud behaviour. We compare our numerical simulations with the experiments performed by Winklhofer et al. [11]. This experimental set up is chosen due to the simple sharp-edged geometry of the entry of the nozzle channel as well as its frequent use in the literature as a benchmark case for computational fluid dynamics (CFD) prediction of cavitating flows. Moreover, although cloud dynamics are not explicitly discussed in this experimental work, it is one of the few experimental data sets that presents qualitative (shadowgraph based vapour fraction) as well as quantitative (pressure and velocity field) measurements. In previous CFD studies of internal geometries relevant to cloud dynamics, the focus was mostly on erosion prediction and geometries like the one from Franc [20,32] were used. However, in the studies of these geometries, limited quantitative comparison is provided with experiments which makes questionable the conclusions drawn only from CFD predictions for the mechanisms of cloud break-off. Moreover, in the investigated cases we used the Linear-upwind stabilised transport (LUST) and MUSCL differencing schemes (the former being a blend of the centred and linear schemes) at pressures relating to the onset of cavitation and the critical cavitation, with different numerical boundary conditions. The introduction of the wave transmissive outlet boundary condition is compared to the typical boundary condition of a fixed pressure value commonly used in previous studies. The physical length of the chamber is also varied to understand better how dependant the flow is, from a numerical point of view, on the geometry.
After having established the extent of the validity, as well as the limitations, of our framework for both cases of cavitation inception and critical cavitation conditions, in the second part of our paper we turn our attention to the mechanisms that govern the cloud break-off dynamics. Using the Winklhofer geometry we have the opportunity to examine the cloud dynamics in conditions close to "choking" that have not been investigated before by numerical studies performed in the same geometry [12,14,15,27]. Previous studies predicted the shape of the free surface of the cloud to be largely flat and generally parallel to the wall for large pressure differences, due to low-order discretisation schemes used, which is inconsistent with experiments. This result however implied a cloud break-off mechanism that was reliant on the re-entrant jet dynamics. In contrast, in our study using a higher order scheme, MUSCL, we predict a vapour cloud topology converging towards the centre of the nozzle which is more consistent with the evidence of the shadowgraphic images. Using these predictions, a detailed analysis of the mechanism of the cloud breakoff is presented. We conclude that the mechanism is dominated by the underlying flow pressure dynamics rather than re-entrant jet dynamics.

Experimental Set up
The overall experimental set-up, along with an enlarged view of the nozzle under investigation, is shown in Figure 1. A micro-channel nozzle etched into a 0.3 mm deep channel was sealed by two transparent windows to allow for the visualisation of cavitation. The nozzle was orthogonal, 1000 × 284 × 300 µm, with the inlet to the nozzle having a radius of 20 µm which corresponds to the "U" geometry described in [11]. The outlet pressure was varied (1.5 MPa, 2.9 MPa, 3.3 MPa, and 4.5 MPa) while the inlet pressure was constant at 10 MPa. The lengths of the exit and entry chambers were approximately 24 mm in the experiment: however, to include the entire geometry within the solution domain, while at the same time resolving the flow in the nozzle, would be computationally demanding. Thus, to model this nozzle accurately, it was desired to shorten the extent of the computational domain representing the lengths of the entry and exit chamber as much as possible. We thus performed a systematic variation of the exit chamber lengths, as well as of the outlet boundary conditions using the LUST divergence scheme and then repeated the investigation with the TVD MUSCL scheme. Table 1 details each case that was simulated. OpenFOAM was chosen as the baseline code [33]. The following section presents the governing equations, along with details of the divergence schemes, phase change modelling, boundary conditions and turbulence model. Figure 2 shows the locations of the various computational probes around the cavitation and cloud break-off regions taken at high sampling frequency (approx. 100 MHz) and spanwise plot location over the width of the nozzle supplemented by Table 2.

Multiphase Modelling
A one-fluid approach was adopted, based on the HEM model [12,34] for phase change and a barotropic closure for the equation of state developed by Kärrholm [12] called cavitatingFoam. Liquid and vapour are treated as one continuum, so that each phase does not require independent transport equations to be solved. Both phases are considered compressible, with the temperature assumed to be constant. The one-fluid approach is a reasonable approximation here due to the high pressure differences obtaining in the experiment. The disadvantage is that the effect of interfaces is not included in the original model formulation. Our previous research [13] has shown that this might be a reasonable assumption for the Diesel fuel that is used in the current work and can be problematic only in cases of fuels/fluids with low viscosity, such as kerosene. A flow chart detailing the solution procedure can be found in Appendix A.

Governing Equations and Compressibility
The governing equation for conservation of momentum in the LES context is as follows Although the experiment was run under quasi steady-state conditions, all flow features and, in particular, cavitation structures, were expected to be unsteady with short time scales. LES is based on the idea of computing the large, energy-containing eddy structures (filtered quantities) which are resolved on the computational grid, whereas the smaller, more isotropic, sub grid structures (sgs) are modelled. In the above equation (and the following ones) the over bars and tildes represent the spatially filtered and density weighted filtered values with the filter width respectively. The filter width is taken as the cube root of the local grid cell volume. Within the current formulation, the only term that maintains the explicit coupling of the effect of small scales to the large scales is τ sgs ij which is determined by the sub-grid scale WALE turbulence model [35], with the LES-dependent constant for the eddy viscosity term set to 0.158. This model has been shown to perform well under similar conditions for Diesel injection [36].
Since a one-fluid approach was used, in order to determine whether the fluid in a cell is liquid, vapour or a mixture, a phase fraction indicator function was solved, namely the vapour fraction: where ρ v,sat = ψ v p sat and ψ v is the compressibility of the vapour. The model assumes that homogeneous mixing between the two phases occurs instantly and thus is better suited for the simulation of flows where cavitation bubbles disperse rapidly such as turbulent cavitation clouds in high velocity flows. It is not so well suited for stratified or slug flows as the two phases may have highly differing velocities requiring better interfacial treatment where a transport equation based on vapour fraction is needed. A preliminary comparison of the HEM approach with the cavitation mass transport approach for different cavitation numbers has been performed by Vogiatzaki et al. [37] which showed that HEM was better for high-speed flows although for slower flows the results from the two models were very similar. Furthermore, the homogeneous mixing assumption avoids the requirement of reconstruction of the liquid vapour interface, as each computational cell is assumed to be fully mixed.
The equations of motion were closed using constitutive relations for the density ρ and the dynamic viscosity µ with the aid of the previously introduced vapour fraction variable α.
The mixture viscosity was given by: For the calculation of the density, using the barotropic assumption with a linear equation of state leads to: The first term governs the liquid's density when α is low. The term ψ is the compressibility of the mixture which is defined as ψ = 1/c 2 with c being the speed of sound of the mixture. If the fluid is fully in the vapour phase, the second term becomes dominant. The property ρ 0 l is given by Other alternative models for compressibility are the Chung [38] and the Wallis model [39]. The literature [12] suggests that using the Wallis compressibility model reproduces the speed of sound (compressibility) in a bubbly mixture in a more physical way, but can be numerically unstable when used for high speed flows. The Wallis model does not assume uniform velocities across the two phases and thus may be better suited for stratified or incompressible mixtures. The Chung method is derived from source terms which include surface tension and bulk modulus and has been shown to work well for dispersed bubbly flows in the region α < 0.3. Figure 3 shows a comparison of speed of sound that can be calculated using the three compressibility models and it should be noted that both the Wallis and Chung models calculate local speed of sound to be, correctly, considerably lower for a mixture than for either the vapour phase or liquid phase. The linear model provides a relatively high speed of sound of the bubble mixture. Although the calculated value of the speed of sound is not used directly in our numerical framework it is important to underline this limitation since this might affect the predictions of areas with Mach number higher than one, as we will demonstrate in the following analysis (see Section 3.2).

Boundary Conditions
The outlet boundary conditions here are considered on the numerical basis of how they will influence the flow, more specifically the cloud break-off mechanism, within the domain. Travelling pressure waves from the violent collapse of the vapour cavities propagate throughout the domain where the large gradient in fluid properties will interact with the outlet boundary. Two boundary conditions were selected to produce opposing scenarios; one where the boundary condition attempts to replicate the pressure gradient so that the waves are not reflected back (referred in the remainder of the paper as "wave transmissive" condition) and the other in which pressure waves towards the boundary are reflected back into the domain (fixed pressure). A fixed value pressure condition was used to force reflection of pressure waves at the boundary. This is not the most optimum pressure condition for pressure driven flows, however is used to help us highlight the scenario that artificial waves interfere with the case physics. In order to avoid the interference of artificial waves that are created in cases that the numerical domain is shorter than the physical domain, the acoustic characteristics of the region beyond the numerical domain needs to be included. The wave transmissive boundary condition is based on a simplified form of Poinsot & Lele non-reflection multidimensional boundary condition called the Navier-Stokes Characteristic Boundary Condition (NSCBC). This boundary condition attempts to identify amplitude variations based on a distance downstream from the boundary to consistently set the boundary value to a given incoming pressure wave [40]. This boundary condition is partially reflective, depending on the distance from the real far field boundary. The shorter this distance is set, the more reflective the boundary is. This boundary condition allows for geometries such as ducts and pipes to reduce cell population at the outlet, as the distance between the boundary and the main flow can be shortened. For this case the far-field distance was set based on the difference between the numerical boundary to the physical boundary from experiment.

Differencing Schemes
The effect of the numerical differencing of the divergence terms in the solved equations is of interest in this investigation. A strategy to achieve results with good accuracy whilst suppressing numerical artifacts is to use a blended scheme with a fixed value blend ratio. The LUST divergence scheme is a fixed blend of two differencing schemes, with the blend fixed at 75% upwind and 25% central differencing [41]. The scheme offers an accurate and stable solution that is computationally inexpensive, however it is not fully bounded. The higher order central differencing contribution is the source for the unbounded solution. Thus, it would be advantageous to dynamically suppress this contribution at regions with large gradients. Total Variable Diminishing (TVD) schemes can offer this option since they represent a class of switching face limited schemes that provide a blended solution with a dynamic limiter. This can suppress the high order contribution at regions of large gradients (reducing numerical instabilities) but allows for the high order solution at smoother regions (increasing the accuracy). As long as the TVD criterion is satisfied, the limiter ensures for a stable and bounded solution. The scheme that we will evaluate in this work is called MUSCL and it is a TVD scheme that uses the monotonised central face flux limiter [42]. The calculation of a dynamic limiter as in the case of MUSCL is computationally more expensive than using a fixed limiter as in the case of LUST, however allows for some high order differencing even at regions of high gradients and can resolve discontinuous functions well, both aspects important for cavitating flows.

Thermodynamic Properties
The solution procedure is described in Appendix A and Table 3 shows the thermodynamic properties that have been used. Table 3. Thermodynamic properties of Diesel [43]. The selection of these properties was based on the review carried out in our previous work [14]. The choice of the individual fuel properties compared to the real physical properties of the diesel fuel have to be considered. The difficulty in the use of Diesel is that it is a blend of hydrocarbons that can vary seasonally as well as geographically and hence has a wide tolerance on its permissible properties. For example in the UK, the specification of common pump diesel must conform to EN590 [43] (which has detailed requirements around sulphur content and bio-fuel components due to their effect on emissions) so that the stated density range is 820 to 845 kg/m 3 . Most authors therefore also choose a value within this range. The values for (ρ v ),(ψ l ) and (ψ v ) however are more difficult to determine since they are not stated in the fuel specifications and include assumptions around dissolved gas content and the exact vapour composition. In this case, values for (ρ v ) range from 0.016 to 6.5 kg/m 3 and 0.1361 was chosen as this is derived from the compressibility of the vapour at room temperature. The liquid vapour compressibility in turn is an estimation based on the ideal gas law as in Kärrholm et al. [12].

Mesh
The geometry was modelled using a structured mesh of 1.1 M elements with a grid density in the nozzle of 506 elements/mm and a depth of 0.3 mm. The minimum cell thickness within the nozzle was 4 µm. The Kolmogorov length scale was estimated to be approximately 0.47 µm and the Taylor scale was of the order 1 µm. Although both of these scales are more meaningful in the context of single phase flows, still they provide an indication of the accuracy of our resolution for the two phases examined in this work. The dimensionless wall distance y + was calculated to be 2.4 so turbulence was handled implicitly rather than requiring explicit wall functions. [14]. Figure 4 shows the numerical grid used for sensitivity analysis and, as can be seen, while the domain length was varied, the mesh quality in the upstream or nozzle domains was kept the same.

Results and Discussion
In this section we initially discuss the sensitivity of the flow dynamics to the differencing schemes (i.e., LUST and MUSCL) of the divergence terms and results are presented for two different outlet pressures, one representing the onset of cavitation and one representing the critical cavitation. We focus on exploring (a) whether the suggested schemes preserve the large pressure gradients and (b) how the prediction of pressure dynamics affect the cavitation development. Alongside experimental comparisons, a discussion on the flow sensitivity to the exit boundary conditions is included in an effort to separate "real" pressure disturbances from "artificial" waves associated with the imposed boundaries and the domain size (i.e., waves reflected from computational boundaries). A link between the numerically reproduced local pressure dynamics and the cavitation cloud break-off mechanisms is provided. Cloud frequencies, as well as vapour fraction fluctuation frequencies, are presented and compared. We also include a short discussion around the choked flow characteristics that some of the operating conditions examined exhibit.

Differencing Scheme and Boundary Condition Sensitivity Analysis
The average mass flow rates are presented in Figure 5 for all the cases that are investigated. The horizontal black line indicates the experimental data. It should be pointed out that at p di f f = 7 MPa and above, the mass flow rate stabilises indicating a "choked" flow. This explains the single experimental value used for comparison for all cases. The numerical predictions are, overall, in fairly good agreement with the experiments especially for the 7.5 MPa case. Under-prediction is mostly noticed for the LUST scheme at 8.5 MPa which is associated with the under-prediction of the vapour fraction at the centreline of the nozzle as we will see in the next paragraphs. It is interesting to note that for p di f f = 8.5 MPa, the inclusion of the wave transmissive boundary condition has a large effect for the case of LUST scheme, while the MUSCL scheme is rather insensitive to the boundary conditions and the size of the geometry.   Table 1) together with the experimental profile. It can be seen that all cases at this pressure condition fit the experimental profile well although all cases also display a slight under-prediction of the pressure around the mid-length (displacement of 1.5 mm) of the nozzle. Table 4 shows the Root Mean Square Error against the experimental data for injection pressures of 7 MPa. The cases with the wave transmissive boundary condition-regardless of the differencing scheme-show a better recovery but this occurs further downstream of the nozzle. At the selected saturation pressure (p sat ) of 0.0054 MPa, no phase change occurs at the nozzle centreline and thus this comparison mostly reflects the ability of the differencing schemes and boundary condition set up proposed to reproduce the underlying flow dynamics. Figure 6. The data correspond to p di f f = 7 MPa, 12 mm computational exit chamber length. Left: Contour plots of the simulated vapour fraction, as well as shadowgraph pictures from the experiment [11]. Right: Pressure profiles along the centreline of the nozzle, including experiment [11]. The notation wT and fP corresponds to wave transmissive and fixed pressure boundary conditions respectively. The two vertical black lines denote the nozzle entry and exit respectively.  Figure 6 (left) shows experimental and computational results for the vapour fraction contour. The results correspond to the wave transmissive boundary condition. The results for the fixed pressure boundary condition cases look very similar and are thus not included here. Cavitation is present at this condition but it is limited to the wall close to the channel entrance and critical cavitation has yet to be reached. The cloud break-off observed at this condition is limited and has minimal influence on the overall flow dynamics. Both schemes under-predict, the stream-wise cavitation extent and the magnitude of the void fraction although the MUSCL scheme appears better.
Critical cavitation based on experimental observation is reached when p di f f exceeds 8.5 MPa. Figure 7 on the right shows the calculated pressure profiles along the centreline of the nozzle for p di f f = 8.5 MPa cases for the 12 mm long computational downstream chamber. Experimental data for this pressure difference are not available for the centreline pressure distribution as they were for Figure 6. One more chamber length (5 mm) is also included since, as shown in our previous work [14], shorter geometries for critical cavitation regime make the flow inside the nozzle more sensitive to the selection of the downstream boundary location and hence boundary conditions. The effect was found to be less pronounced at lower pressure differences (for example at p di f f = 7 MPa) when cavitation is no longer critical. Using these two cases we are able to analyse better the dependency of the pressure dynamics on the domain geometry and the outlet boundary conditions. The majority of the cases compared predicts similar pressure profiles apart from the 5 mm geometry using the LUST scheme with fixed pressure boundary conditions which shows an increased pressure.  [11]. Right: Pressure profiles along the centreline of the nozzle. The notation wT and fP corresponds to the wave transmissive and fixed pressure boundary conditions respectively. The two vertical black lines denote the nozzle entry and exit respectively.
A corresponding colour contour plot of the vapour fraction is presented in Figure 7 (left) for wave transmissive outlet conditions that based on Figure 7 on the right it is the boundary condition that makes the results more insensitive to the geometry size as expected. At critical cavitation conditions, the experimental data shows that the vapour fraction topology reaches almost the centre line of the nozzle where the flow becomes "choked" (following the experimental paper term) [11]. It should be noted that in reality it is difficult to fully determine the nature of the cavitation cloud at the nozzle centre from the shadowgraph images alone since these are line of sight and represent the ensemble average beam intensity from 20 images as a 'probability' (no calibration was carried out to determine the vapour fraction required to reduce beam intensity to zero). In this context, the indications from the shadowgraph data are corroborated by mass flow measurements showing that the flow is indeed choked because of the stabilisation of the exit mass flow rate. The LUST scheme shows a cloud topology with the cavitation region appearing smooth and following a vena contracta-like curvature. In comparison, the MUSCL scheme shows a more turbulent vapour region and, possibly, cloud break-off features. As the MUSCL scheme is better able to resolve large gradients, it is able to capture local inhomogeneities whilst the LUST gives some artificial smoothing across the interface between the wall-attached vapour cavity and the fluid moving along the centreline. Qualitatively, this model captures some part of the experimentally observed extension of the cavitation cloud towards the duct's centreline. Figure 8 presents 3-dimensional snapshots of the instantaneous liquid vapour interface (for an iso-surface of vapour fraction of 0.5) predicted by the MUSCL and LUST schemes at p di f f = 8.5 MPa using the 12 mm long computational chamber geometry with the wave transmissive pressure condition. The figure is included in order to provide more insight into the differences in the predictions of the cavitation cloud topology when different schemes are used. The LUST scheme predicts a re-entrant jet whereas this is not evident in the MUSCL scheme which shows much greater wrinkling of the interface as it interacts with the flow vorticity. For completeness we have included also the plots for p di f f = 7 MPa in order to demonstrate the evolution of the cavitation area when the pressure difference is increased. The stream-wise pressure distribution and the vapour fraction along the near wall region for all cases with the fixed pressure and the wave transmissive condition are presented in Figure 9. The plots on the top corresponds to the stream-wise pressure distribution along the near wall region for all cases alongside experimental data at 7 MPa. The plot on the bottom correspond to the vapour fraction profiles close to the wall calculated along the same line as for the pressure profiles. Although experimental data are not available, the vapour fraction figures are included in order to demonstrate the link between the predicted pressure and the vapour cloud. For the p di f f = 7 MPa case, the experimental data indicate a sharp decrease in pressure at the channel entrance-essentially a Bernoulli effect-but also a very quick recovery which indicates that the cavitation cloud has either already detached or that bubbles will start quickly collapsing. Both MUSCL and LUST differencing schemes capture this behaviour reasonably well for p di f f = 7 MPa and the 12mm long computational chamber. The recovery shown is slightly under-predicted compared with experiments, but better recovery is achieved compared to previous numerical work [44]. Table 5 shows the RMSE against experimental data. The predictions are better for the wave transmissive condition, and as seen with the centreline, the smallest error was found when coupled with the MUSCL scheme. At p di f f = 8.5 MPa the MUSCL scheme appears to be rather insensitive to the downstream length regardless of the outlet boundary condition (observation consistent with the previous figures) while the LUST scheme is more sensitive when fixed pressure outlet boundary condition is used. For the vapour profiles, at p di f f = 8.5 MPa, it can be seen that there is strong cavitation close to the wall over more than half the upstream length of the nozzle. Part of this is 'sheet' cavitation (values close to 1) while the rest is where the cloud is starting to form. In the downstream half of the nozzle for the fixed boundary condition, no cavitation is present close to the wall, implying that the cloud has already detached and potentially collapsed. When the wave transmissive outlet is used, the LUST scheme (regardless of the length of the domain) predicts some vapour structures at the downstream half of the nozzle. This is consistent with what was observed in Figure 7 (left) that indicate the inability of the scheme to correctly reproduce the vapour detachment from the wall and convergence towards the nozzle center.  In order to highlight the interplay of the vapour cloud with the underlying flow field dynamics Figures 10 and 11 are included. Figure 10 shows the colour-coded, instantaneous distributions of the flow vorticity (taken from the velocity of all components) for LUST and MUSCL in the symmetry plane of the nozzle (12 mm long computational chamber). For LUST cases, the vorticity is more pronounced along the liquid vapour boundary and the trailing region. Within the area where a re-entrant jet exists, the greater vorticity is at the head of the jet. For MUSCL cases, the vorticity is less concentrated along the liquid vapour boundary. The free stream at the entrance of the nozzle becomes chaotic as it starts to interact with the nozzle boundary and the cavitation interface. The greatest extent of the vorticity can be observed along the upper and lower boundary of the nozzle both within, and out, of the vapour cavity region. This hinders the development of the re-entrant jet and any upstream flow due to the vortical structures present. This conclusion is also consistent to what can be observed in Figure 8, although for the LUST scheme the re-entrant jet was present while for the MUSCLE scheme this was not so obvious.    Figure 10. A 0.5 iso contour of the vapour fraction is superimposed. It can be seen that for the LUST scheme, the flow at the nozzle centre is laminar while for the MUSCLE scheme turbulent structures are present in the centre. Thus, it could be argued that the mechanism for cloud break-off in conditions close to choking seems to be associated with turbulence and subsequent wrinkling of the surface which eventually leads to break-off in a way similar to how liquid atomisation occurs, rather than the re-entrant jet that is suggested in the literature. It has also been reported in the literature [45] that nozzle turbulence is linked to the collapse of the cavities and that a small change in the size of the cavity can cause a substantial increase in the turbulence level and momentum thickness in the boundary layer downstream. For the LUST scheme, since the cavities always remain very close to the boundary layer, the turbulence generated by the collapse is dissipated quickly. For the MUSCL scheme, the cavitation extends towards the centre of the nozzle and the vorticity generated by the collapse is maintained by the high speed core fluid (where the velocity is above 100 m/s). Figure 12 presents streamlines projected along the centre of the nozzle for LUST and MUSCL. Recirculation, due to stagnant flow from the presence of the re-entrant jet, can be clearly seen in the LUST case with the flow stagnating around this region (blue lines). For the MUSCL case, the turbulent structures prevent the re-entrant jet from forming coherently, where the flow starts to curl and oscillate downstream from the cavitation cloud break-off. Reverse flow (negative velocity) is observed in small pockets along the near wall but are not large enough for any jet structures to form.

Wave Visualisation and Flow Choking
The numerical modelling of cavitation within a compressible framework (mostly density-based schemes) can, in theory, capture all the waves generated in the computational domain by the flow characteristics and the phase change phenomena. However, in practice it is hard to separate the "waves", or disturbances, that appear due to numerical reasons (such as reflecting boundary conditions and/or artificial pressure oscillations) from the real pressure waves within the domain from, e.g., bubble growth and collapse. Moreover, another characteristic of flows that involve phenomena dependent of pressure waves that is desirable to reproduce numerically is that the energy of these pressure waves dissipates relatively quickly with distance. However, the dissipation rate in simulations is controlled usually by the available numerical schemes. A wave can be artificially smoothed out by low-order numerical effects (due to numerical dissipation), or spurious oscillations near the wave surface might arise if a high-order numerical scheme is used (due to Gibbs phenomena). Using a higher order scheme such as MUSCL could therefore augment the propagation of these waves if incorrect boundary conditions are used. A third challenge is the difficulty in modelling the speed of sound in a multiphase fluid where the size, composition, distribution and shape of cavitation bubbles all play a part in determining the value of the local speed of sound. In the current work, for example, we have used a simple linear relationship for the compressibility (and consequently the speed of sound) due to stability considerations rather than accuracy. This results in a higher speed of sound in the mixture part of the flow, leading to predictions of local Mach numbers lower than expected. In this section we will provide more insight into the link between the boundary conditions and numerical schemes with the pressure oscillations that govern the mechanisms of the cloud break off for conditions close to choking. Figure 13 shows the instantaneous pressure distribution in the region immediately downstream of the nozzle for the 12 mm computational chamber length with p di f f = 8.5 MPa. Cases for both the fixed pressure and wave transmissive boundaries are represented at three time steps 0.05 ms apart, chosen to highlight the interaction between the propagating waves and the outlet boundary. A high pressure wave that was generated by a disturbance in the centre of the free-stream can be seen in the case that a fixed pressure outlet was used. This wave is reflected by the upper and lower boundaries and is propagating towards the outlet boundary (upper frame). The middle frame shows that this wave is reflected back into the domain, however it has now been inverted (i.e., it was a compression wave and is now a rarefaction wave). The lower frame shows how this influences the dynamics of the domain as it interacts with the high pressure oscillations around it. This is the expected physical behaviour if the outlet boundary was an open pipe into a reservoir. For the case with the wave transmissive outlet boundary condition, similar disturbances central to the free stream flow were harder to find. The scale used for all the figures, related to the wave transmissive outlet conditions, demonstrates the overall lower pressures seen in this region using this boundary condition though it is difficult to directly attribute the drop in pressure to the suppression of reflected pressure waves at the outlet. A compression wave originating in the upper part of the domain near the outlet boundary can be observed in the upper frame. As the wave propagates away from its point of origin, it is reflected back by the solid upper boundary but is transmitted away (rather than reflected back) by the outlet boundary. Thus does not affect the rest of the flow closer to the nozzle as strongly as in the case of fixed boundary conditions.
A similar representation of the pressure field for the nozzle and immediate nozzle exit region can be seen in Figure 14. The same pressure scale was used for a comparative analysis around this region with an additional blue contour which indicates a vapour fraction of 0.5. The superposition of reflected waves creates an area of higher pressure which, for the fixed pressure case, coincides with the centreline pressure distributions as seen in Figure 7. The initial cavitation collapse occurs within the nozzle, with waves being generated inside the nozzle region. The overall cavitation around the exit of the nozzle is suppressed, as observation of Figure 9 confirms. Conversely, with the wave transmissive boundary condition, this collapse occurs outside the nozzle exit region. In the remainder of the section we will inquire more closely as to what is the dominant cloud break off mechanism and if it is consistent with the two mechanisms currently identified in the literature (bubbly shock or re-entrant jet). The bubbly shock mechanism requires that the local speed of sound reaches a low enough value, and thus a high enough Mach number, for it to be the dominant mechanism for cloud break-off [5]. As travelling pressure waves from collapsed cavities pass through a vapour cloud, the speed of the wave is much closer to-if not greater than-the speed of sound within the mixture, causing the shock to propagate inside the cloud. Figure 15 is a colour contour plot of the Mach number across the nozzle geometry which was calculated from the speed of sound based on the linear mixture compressibility. For both the fixed pressure and wave transmissive cases, high Mach numbers can be seen within the cavitation cloud extent, with a concentration along the α v = 0.5 boundary, but the value never exceeds 0.26. In order to understand the extent of the discrepancies of our prediction of the speed of sound we use as reference a different approximation reported by Shamsborhan et al. [46]. The speed of sound within the cavitation region can be estimated for P (di f f ) = 8.5 MPa, the point at which critical cavitation was reached experimentally, to be c mixture ≈ 31 ms −1 . This is much lower than the reported approximate free stream velocity of 120 ms −1 which would lead to an estimated Mach number within the cavity region to be Ma ≈ 3.9. According to Ganesh et al. [5], for bubbly shock to become the main mechanism, the Mach number should exceed unity and a value of ≈ 3.1 was found to be the critical value. Thus, if a better speed of sound approximation was used then the calculated local Mach numbers of Ma ≈ 3.9 would be consistent with the bubbly shock formation. Another aspect that limits the minimum local speed of sound within a bubbly mixture is that the pressure cannot be below the vapour saturation point, P v,sat = 5400 kPa based on the HEM model. This means that, at this point, the flow is no longer fully compressible at constant temperature, and could further decrease the local speed of sound though this is difficult to calculate. Looking at the cavitation numbers, the critical cavitation number for the point of dominant bubbly shock for the hydrofoil case of Ganesh et al. [5], was approximately 1.95, which is higher than the value of approximately 0.65 reported by Jahangir et al. using internal venturi geometry. The cavitation number in our case, from Table 6, is 0.16 which also indicates a cavitation flow of high intensity. Although in our study the bubbly shock traversing a mixture is not fully captured because of the over prediction of the mixture speed of sound, the lack of re-entrant jet in the experimental evidence was reproduced accurately by our simulations. Summing up, this suggest that for p di f f = 8.5 MPa it is reasonable to expect that the main mechanism for cavitation cloud break-off to be bubbly shock and this mechanism is different to the re-entrant jet mechanism demonstrated in previous studies for the same case.

Bubble Cloud Dynamics
The importance of identifying the dominant frequencies of the phenomenon of the periodic break-off of cloud cavitation from attached wall cavitation is linked to various applications of interest among them the erosion potential of cavitation. To analyse the break-off frequencies in our study, seven "numerical" probes were placed at different locations around the liquid-vapour interface to measure the fluctuating pressure. Figure 2 shows the various locations of the probes in the near wall region of the nozzle and Table 2 gives the corresponding coordinates. A Hanning window was applied to the input signal before an FFT was performed to analyse the frequency response at the probe locations. The cloud cavitation break-off will not be the only aspect of the flow that will be generating a fundamental frequency. The geometry of the nozzle and the chamber may produce a resonance (perhaps several) frequency that is dependent on the topology. These waves are reflected back by any solid boundary and will interfere with other pressure waves present, as also reported by [19]. In their work (albeit using a different geometry) the higher frequencies among the four reported (1 × 10 3 Hz) were attributed to this phenomenon.
To understand better the geometric influence on the frequencies in this study, a simplified two dimensional test case was also performed. Figure 16 presents the setup for this test case, comprising of three walls and an outlet boundary that were set using the same conditions discussed earlier. The initial pressure field was distributed with a high and low pressure regions both sides of a pressure wave front. The cell density of the two dimensional mesh corresponds to the same density as that of the nozzle geometry in order to reduce the dependence on the grid.  In order to be able to analyse the frequencies (and to resolve to the lower frequencies reported in literature [6,21,47]) the computational time at a quasi-steady state needs to be as long as possible: however this increases, considerably, the computational cost. Thus, a single case, (MUSCL with wave transmissive outlet boundary condition with 12 mm long computational downstream chamber at p di f f = 8.5 MPa), was selected to run for an extended period of 5 ms resulting in the spectral plots of the pressure and α signals presented in Figure 17. This extended case has a minimum resolvable frequency of 600 Hz. Both pressure and vapour cloud frequencies are reported. A dominant frequency in pressure distribution exists around 2000 Hz, with harmonics at 4000 Hz and 8000 Hz. The reported frequencies for the two dimensional geometric test presented in Figure 17 suggest that the resonant frequency due to the geometry corresponds to approximately 1500 Hz and a harmonic of approximately 2700 Hz. These frequencies can be identified within spectrum from the cavitation with a lesser intensity to the dominant frequency.
This shows that the source of the dominant frequency is not geometric but rather stems from the flow dynamics. These frequencies can be identified in the 12 mm 8.5 MPa case, with lower amplitude than the flow induced frequencies. For the vapour fraction α, the dominant frequency peak is also at 8000 Hz with the other harmonic frequencies at a little lower amplitude. This results in a Strouhal number of approximately 0.053 which is slightly lower than expected for cavitating flows but within the range seen in relevant literature [6]. For the vapour fraction, the magnitude of the higher frequencies tapers off before increasing gradually with a peak at approximately 300,000 Hz. The high temporal sampling frequency of the probes allow for the identification of these very high frequencies due to the collapse of very small cavities, however these frequencies are not captured in the pressure field. This shows that although the cavitation break-off frequency and the pressure oscillations have a strong link are not entirely related. The Kolmogorov turbulent timescale was estimated to be τ ν ≈ 5.39 µs. Thus, the very high frequencies observed with volume fraction distribution above ≈200,000 Hz, are not translated to the pressure field and can not be directly attributed to the generated turbulence. We believe these high frequencies are present due to numerical discontinuities that exist at the very small scales for multiphase flows.

Conclusions
The aim of this work was to identify the interaction of pressure dynamics and the cavitation bubble cloud break off mechanisms in a micro channel for various pressure differences including conditions approaching flow choking. Two different differencing schemes (a second order TVD scheme called MUSCL and the LUST scheme) were compared in terms of their ability to maintain the steep flow gradients. The predicted pressures, velocities and vapour distributions at cavitation inception and the critical cavitation operating points were compared to experimental data. The exit chamber downstream geometry was varied from 5 mm to 12 mm, with and without a wave transmissive outlet boundary condition to also asses to what extent the pressure waves present in the domain represent physical waves or numerically dependent oscillations. The main conclusions of our work can be summarised as follows:

1.
Cloud break-off for LUST cases is driven by the presence of the re-entrant jet. This cloud break-off characteristic is expected in cavitating flows with low intensity. MUSCL displays characteristics of cloud break-off driven by turbulent wave propagation from the collapse downstream of cavities with the re-entrant jet effectively being overcome by these effects.

2.
The introduction of the MUSCL scheme produced a more accurate cavitation cloud topology, when compared to experimental data, which converged towards the nozzle centreline at choked conditions. This was observed consistently for both the 5 mm exit chamber and 12 mm exit chamber geometries showing that the produced cavitation becomes largely independent of the downstream geometry. The same was not true for the LUST cases where a noticeable change in the centreline pressure distribution was observed changing from 12 mm to 5 mm exit chamber length.

3.
We demonstrated that interference caused by reflective pressure waves from the boundaries affects the dynamics of the cavitation cloud break-off. Higher pressures were observed around the nozzle exit for fixed pressure condition as expected.

4.
The wave transmissive boundary condition had a considerable effect on the flow. It reduced the majority of the reflected shock waves and has a lower mean pressure in the region the outlet boundary and the nozzle outlet compared to the fixed boundary condition. The use of the wave transmissive condition also made the results largely independent of the geometry size.

5.
Although some of the characteristics of choked flow described in the experiments appear to be reproduced by the model (vapour extent, radial velocities and mass flow rates), the linear compressibility model used prevents other physical aspects of fully choked flow such as local speed of sound to be accurately captured.

6.
A frequency analysis was undertaken and compared with other similar works indicating that cloud shedding and bubble collapse frequencies are present. However, due to the spectral range and thus computation expense required, this could only be done on one case for lower spectral bounds.
Overall it has been shown that the dynamics of cavitation are not linear due to the complex interfacial-boundary interaction and the compressible characteristics by the collapse of vapour cavities. The propagation of pressure waves were better realised with the MUSCL differencing scheme and wave transmissive outlet boundary condition. This led to the re-entrant jet being overcome by a bubbly pressure wave mechanism. With these considerations made, a frequency near that expected from cavitation collapse as well as other generated harmonic frequencies were observed.