Compressible bubble dynamic simulations with central-upwind schemes

This paper discusses the implementation of an explicit density-based solver, based on the central-upwind schemes originally suggested by Kurganov, for the simulation of cavitating bubble dynamic flows. Explicit density based solvers are suited for highly dynamic, violent flows, involving large density ratios, as is rather common in cavitating flows. Moreover, the central-upwind schemes have the advantage of avoiding direct evaluation of the Jacobian matrix or estimation of the wave pattern emerging from Euler equations. Second order accuracy can be achieved with TVD MUSCL schemes. Basic comparison with the predicted wave pattern of the central-upwind schemes is performed with the exact solution of the Riemann problem showing an excellent agreement. Then several different bubble configurations were tested, similar to the work of Lauer et al. (2012). The central-upwind schemes prove to be able to handle the large pressure and density ratios appearing in cavitating flows, giving similar predictions in the evolution of the bubble shape.


Introduction
The subject of cavitation bubble dynamics has been studied traditionally using the Rayleigh-Plesset equation [1,2]. However, either the original Rayleigh-Plesset equation, or its more complex variants, assume that the bubble shape is perfectly spherical. In practice this is not the case, since many works, both experimental [3] and numerical [4][5][6], suggest that the bubble shape may be strongly deformed in the presence of pressure fields (e.g. due to gravity, due to passing sound waves), or due to the presence of boundaries (walls, free surfaces, etc.). Thus, if one wishes to predict the asymmetric bubble collapse near a wall, then it is necessary to do so by properly integrating the Navier Stokes equations. The complexity of the involved flow pattern is significant, since the flow is multiphase, involving a strongly deforming liquid/gas interface and very high velocities. In the past, there have been efforts to perform such simulations; one of the first was the pioneering work of Plesset and Chapman [4]. In the present investigation, a method similar to the one used by Adams et al. [7] and Pohl et al. [8] is employed; the cavitation bubble is described as a density difference of a single fluid, governed by a complex equation of state which represents the isentropic phase change due to cavitation. As it will be explained in the next section, a central-upwind scheme is applied, since it has some attractive features. The aim is to predict the outcome of several bubble collapse events near wall configurations, in order to determine the performance of the scheme employed.

Numerical method
The Euler equations are resolved with an in-house code in conservation form, considering the influence of cylindrical symmetry, see Toro [9], to reduce the computational cost. A piecewise barotropic equation of state is used, which is a combination of the Tait equation of state above saturation and a formula that resembles the isentropic transformation within the saturation dome [10]: The values used for the equation of state are: density of saturated liquid ρ sat,L =998.2kg/m 3 , speed of sound of the saturated liquid C sat,L =1450m/s, saturation pressure p sat =2339Pa, bulk modulus B=294MPa and the exponent n=7.15, commonly used for weakly compressible liquids. In order to evaluate the flux at the interface of the finite volumes, the central-upwind scheme of Kurganov is used [11]. Linear interpolations are utilized at the cell interfaces, handled with the MUSCL scheme with a SuperBee limiter [9]. This scheme has the advantage of being universal, in the sense that it does not need the tuning of the AUSM+up scheme coefficients, while it does not require the identification of the wave structure of the solution of the Riemann problem, as e.g. the HLLC or Roe solvers do. On the other hand, the 2nd (or higher) order extension ensures low numerical diffusion. Boundaries are handled either as transmissive or as slip walls, depending on the configuration (see Toro [9]). Time integration is performed in an explicit manner, with a splitting scheme [9], i.e. initially for the homogenous part of eq. 1 and then for the source term. In this preliminary work, 1 st order Euler integration is used, with a CFL of 0.5, whereas in the future higher order Runge Kutta integration will be implemented.

Validation of the numerical method
In Fig. 1a the solution of the Kurganov scheme for the Riemann problem with initial conditions ρ L =1002.88kg/m 3 and ρ R =9.99kg/m 3 , u=0m/s everywhere is shown; also the Lax-Friedrichs and exact solution are shown for reference. It is of interest that the central-upwind scheme is successfully capturing the correct wave pattern, with the same spatial resolution as the Lax-Friedrichs scheme without smearing of the shock, due to the inherent numerical dissipation of the latter. It should be highlighted here that obtaining the exact solution of the Riemann problem is not trivial for an arbitrary EOS, such as the one in eq. 1, due to the nature of the Riemann invariants in the rarefaction zones. In order to validate the 2D axis-symmetric solver, which will be used later, a comparison with the 1D solver with spherical symmetry was performed. The initial conditions for this comparison are similar to the above, ρ=9.99kg/m 3 for R<1m and ρ=1002.88kg/m 3 for R>1m, resembling an implosion configuration. In Fig. 1b

Results
The case of interest is the collapse of a water vapor bubble in the vicinity of a wall, in the same arrangement as used in Lauer et al. work [5], using the framework analyzed in section 2. The bubble has a radius of 400µm and its center is placed 416, 140 and -140 µm from a wall. The surrounding fluid has a pressure of 100bar, whereas the pressure within the bubble is equal to the saturation i.e. 2340Pa. In Fig. 2 indicative instances of the bubble deformation during the collapse are shown for the three configurations. The influence of the wall distance is apparent: when the initial bubble centre is above the wall, then the jet effect is formed on the axis of symmetry, causing the bubble to take a heart like shape or a torus in the cases of 416µm and 140µm respectively. On the other hand, in the case where the initial bubble centre is positioned below the wall, the jet forms in the x direction at the wall, deforming the bubble in a pin-like shape. In all cases, the formation of the jet, either on the symmetry axis, or in the x-direction causes the development of very high pressures due to momentum focusing, see Fig. 2. At the late stages of the bubble collapse, the jet will eventually impact on the wall, causing pressures of the order of 10000bar; such pressures are well beyond the yield stress of many common materials (e.g. SS316L has a yield stress of the order of 2-4 . 10 3 bar), implying that such bubble collapse configurations will contribute to the erosion damage of the underlying solid material. From a numerical point of view, the employed scheme performed well, in the sense that it is able to handle pressure ratios of almost 5 . 10 5 and density ratios of 1000, without serious problems. High accuracy enabled a clear capturing of the interface within 1-2 cells, without oscillations, thanks to the TVD properties of the MUSCL scheme.

Conclusion
This paper outlines the development of an explicit density based solver for cavitating flows, based on the central upwind schemes of Kurganov et al. and the application on bubble collapse using 2D axissymmetry conditions. The schemes have been tested in comparison with the exact solver of the Riemann problem, showing excellent accuracy and robustness. Application of the schemes on the bubble collapse cases showed a similar collapse pattern with the one that has been reported by Laeur et al. and similar pressure levels on the wall, even though a Homogenous Equilibrium assumption is used for the thermodynamic model. Overall, the performance of the central-upwind scheme was very good, considering the violent nature of the collapse and the involved pressure and density ratios. One of the main targets in the future is the implementation of higher accuracy in the time marching and possibly inclusion of thermal effects, with a potential application the simulation of bubble clusters.