Numerical study of real gas effects during bubble collapse using a disequilibrium multiphase model

Graphical abstract


Introduction
Thermal effects occurring during the collapse of gaseous bubbles [1], including sonoluminescence [2], air dissociation and chemical reactions [3] are now well documented. While the bulk liquid temperature does not change significantly compared to the inner bubble content [4][5][6][7], the latter can reach enormous temperatures during the collapse, of the order of thousands of degrees Kelvin as computational studies for both spherical [8,9] and non-spherical bubble collapse cases [10,11] as well as molecular dynamics [12] suggest. A precise determination of the bubble thermodynamics is important in different areas such as in sonochemistry [13,14] [1,3] and ultrasound therapy such as Highintensity Focused Ultrasound (HIFU) to ensure safety and efficiency [15] [4]. Numerical models utilised for the simulation of spherical bubble collapse under such extreme condition typically have employed zero-dimensional approximations, such as the Rayleigh-Plesset [16], Keller-Miksis [17], or Gilmore [18] models. Despite their simplicity and low computational cost, their simplifying assumptions such as spherical symmetry and spatial uniformity of the temperature, limit their applicability to describe the flow physics of more realistic configurations. These include, for example, the asymmetric bubbles collapse that occurs in presence of shocks or near solid boundaries. Moreover, homogeneity of the temperature distribution is affected by Peclet number. [19,20].
To account for such effects and overcome the foregoing limitations of zero-dimensional models, multi-dimensional computational methods have been proposed. In this regard, there are different numerical approaches developed for predicting the temporal displacement of the gas-liquid interface, namely interface tracking [21][22][23][24] or interface capturing [25][26][27]. Moreover, some developed approaches [28][29][30] are based on the solution of the zero-dimensional models. Regardless of the numerical approach, the thermodynamic closure utilised in the flow solvers plays a crucial role in predicting the temperature during the collapse. The vast majority of the relevant publications employ the ideal gas EoS (Equation of State); indicative studies include [31][32][33] for spherical and [34][35][36][37][38][39][40][41] for non-spherical collapse cases.
Nevertheless, the ideal gas EoS does not provide accurate estimates of the gas temperature by ignoring ionization and dissociation as well as the dependency of the enthalpy on pressure which affects the compression heating or decompression cooling. There are only a few studies where a non-ideal gas EoS is used to model the thermodynamics of the bubble. Moss et al. [8] developed a 1D solver with spherical symmetry assumption and used an analytical EoS for the air bubble that includes vibrational excitation, dissociation, ionization, and a repulsive intermolecular potential. Extremely high temperatures up to 1.74 × 10 6 K and 1.16 × 10 7 K have been reported with and without considering the air dissociation and ionization respectively. In [9,42] a 1D model of spherical bubble is employed where the vapour bubble thermodynamics was modelled as a hard-core van der Waals gas. It was shown that a 1 mm vapour bubble initially at saturation pressure surrounded by water at atmospheric pressure can reach a collapse temperature above 10 4 K. The only work known to the authors that has considered non-ideal gas EoS in multi-dimensional flow solvers is [10]; a front tracking method was employed for resolving an argon bubble interface motion under the influence of a strong incident shock wave with varying pressure of 0.1-1000 GPa. The authors compared the results with ideal gas and two real gas EoSs namely the quotidian EoS, which utilise the thermodynamic functions from the Helmholtz free energy with the electronic contribution, and the SESAME database [43,44], which considers the formation of plasma and ionization. It was revealed that the real gas EoSs estimate lower collapse temperatures depending on the shock strength; with a temperature difference of 3.5 × 10 4 K and 1.5 × 10 7 K for incident shock pressure of 1 and 1000 GPa, respectively.
In the present work we are using the six-equation method [45] to investigate thermal effects during bubble collapse accounting for the real gas thermodynamics. This method incorporates two distinct equations for the specific internal energy of each phase. Each phase is allowed to have different pressures. Mechanical equilibrium is imposed after the integration of the conservation equations, taking into account the conservation of the mixture's total energy an equilibrium pressure is calculated. Here we present a numerical approach for imposing the mechanical equilibrium that can incorporate real gas EoSs in a tabulated form. Results are presented for three real gas EoS and compared against predictions obtained with the ideal gas EoS. To demonstrate the role of the non-ideal gas effects at different collapse strength, the spherical collapses are simulated at initial pressure ratios in the range of ≈ 7-353, defined as the initial liquid pressure over the gas pressure. Subsequently, a shock-induced non-spherical collapse in an ultrasound field of a lithotripter is simulated where the collapse pressure reaches the order of GPa.
The first aim of the present study is to gain insight into the impact of the real gas thermodynamics modeled by the cubic [46,47] and Helmholtz [48] EoSs on an air bubble not only for a spherical collapse but also for a practical non-spherical case in Biomedical Science. For the latter, the wall pressure discussed in the case of using different gas EoSs. At the

Abbrevations
Mass fraction Z Acoustic impedance (Pa⋅s/m 3 ) z axial coordinate axis same time, the second aim is to develop the disequilibrium multiphase method [45] to be compatible for any arbitrary equation of state through tabulated data. The remainder of this paper is organised as follows. The adopted equations of state are described in Section 2. The multiphase method and the adopted numerical schemes are explained in Section 3. The section presents various numerical results Section 4. Lastly, some concluding remarks are given in Section 5.

Equation of state
In this study, the stiffened gas EoS [49] is employed for the liquid phase. For the gas phase, in addition to the ideal gas EoS, three distinct real gas EoSs are utilised, i.e.: (1) the Helmholtz EoS [48], and the cubic EoSs of: (2) Peng-Robinson (PR) [46] and (3) Redlich-Kwong Peng-Robinson (RKPR) [47]. Although the Helmholtz EoS is based on experimental data, its range of applicability does not cover the extreme thermodynamic state of the collapse studied in the current work. Therefore, herein, the cubic EoSs with the benefit of a wider range are also employed. The equations of state are described in A. To compare the behaviour of the gas EoSs, the deviations of temperature predictions among the equations of state for a simple isentropic gas compression are illustrated in Fig. 1 The compression starts from assumed values of 300 K and 1 bar for the temperature and pressure, respectively. It is clear that as the compression ratio (defined as the final over the initial pressure of 1 bar) increases, the difference between the ideal gas and real gas models (cubic and Helmholtz EoSs) increases. Notably, at a compression ratio of just 300, the difference in temperature prediction compared to that obtained using the ideal gas model is beyond 10%. At a compression ratio of 4 × 10 4 , this difference is higher than 35%. Such compression ratios are commonly found in bubble collapse cases. For example, in the study of a single bubble excited by a lithotripter [50], which is also simulated in the current work, the bubble starts the compression with the initial atmospheric pressure and reaches a collapse pressure more than 4 GPa. In another study [9], a laser induced bubble is initially at the vapour saturation pressure while it is compressed to the peak pressure of ≈ 10 GPa during the collapse. Comparing the three real gas models, it is observed that the RKPR EoS behaves practically identical to the Helmholtz EoS. The PR EoS, however, exhibits a ≈ 5% difference at a compression ratio of 4 × 10 4 compared to the other two. Observing this, the RKPR EoS is preferred over the PR EoS for the non-spherical collapse case in this study where the bubble pressure reaches pressures higher than 4 GPa. It is noted in Fig. 1 that at higher compression ratios, comparison between the high order EoSs is not attempted, as the Helmholtz EoS is not reliable beyond its calibration range, giving unrealistic density variation as well as problematic values for heat capacity and speed of sound. This is a direct consequence of the high order nature of the Helmholtz model, which implies that its monotonicity is not guaranteed beyond the calibration range.
In this study, the stiffened gas and ideal gas EoSs are used in their parametric form due to the simplicity of their mathematical expressions. The real gas models, however, are implemented through tabular form. This offers a more versatile framework applicable to any EoS, i.e., the CFD solver does not need to be modified when using different tabulated EoSs. Moreover, it disengages the computational problem of evaluating the EoS from the CFD computations and circumnavigates the problem of deriving explicit solutions from the EoS with implicit expressions.
Each table of the real gas EoSs is a rectangular structured temperature -pressure grid with fixed intervals of T and log 10 p. The considered range for temperature with 121 cells is [60, 17000] K for the cubic and is [150, 17000] K for the Helmholtz EoSs. The pressure with 375 cells ranges in [2300, 1.1 × 10 10 ] Pa for the cubic and [2300, 4 × 10 9 ] Pa for the Helmholtz EoSs. It will be shown that while the properties range for the cubic EoSs covers all the collapse cases in the present study, the Helmholtz EoS fails in simulation of the highest initial pressure ratio case.

Numerical method
In the present work, the bubble collapse is modelled using a so-called 'six-equation model' of [45] which stems from the Baer-Nunziato (BN) model [51]. In the BN model individual momentum and specific energy equations are considered for each phase. Thus, each phase possesses its own velocity, temperature and pressure, i.e., full disequilibrium. The BN model can be further simplified to reduced models [52,53], that consider single velocity (kinetic equilibrium) and single pressure (mechanical equilibrium). However, these further reduced models present issues regarding volume fraction positivity and speed of sound monotonicity along with derivation difficulties for the Riemann solver when considering both phases according to [45]. To overcome these issues, the six-equation model of [45] was developed in which the phases have the same velocity but different pressures and temperatures. According to its developers [45], the six-equation model is not a physical model but a step model for the reduced models to overcome the mentioned issues of [52]. In this model, the phases will reach the mechanical equilibrium through a relaxation process at infinite rate which will be described in the solution algorithm. The model incorporates the mass and the energy conservation equations for each phase, a single momentum conservation equation for the mixture (considered as one equation in vector form) and the volume fraction transport equation for the first phase. In addition, the mixture total energy conservation equation is also solved as the seventh equation to ensure the total energy conservation. More details about this model can be found in literature, e.g., [54,31,33,55,56].
In the present work, the six-equation model is used in the 1D spherical and 2D axisymmetric coordinates (cylindrical coordinates with azimuthal symmetry) to save computational cost: where: Fig. 1. Comparison of the temperature obtained with ideal and real gas EoSs for different compression ratios in isentropic compression (T 0 = 300 K and p 0 = 1 bar).
where s rlx , s nc , and s g indicate relaxation, non-conservative, and geometric source terms, respectively. Moreover, subscripts 1 and 2 denote the water and air phases, respectively. Also, the following notation is adopted: r, z (coordinate axes), t (time), ρ (density), p (pressure), α (volume fraction), u (r-direction velocity), w (z-direction velocity), e (specific internal energy), E (specific total energy). The interfacial pressure p I is defined as follows which is an estimate in the limit of equal velocities introduced first in [57]: where Z k = ρ k c k denotes the acoustic impedance for phase k with speed of sound c k . The value of 2 and 1 for the coordinates switching parameter β correspond to the 1D spherical the r-direction and 2D axisymmetric coordinates in the (r, z) directions, respectively. As can be seen in Eq. 2, this model assumes an initial disequilibrium pressure meaning that each phase has its own pressure. As pressure waves reach an interface, they propagate to the next phase by a pressure coupling between the phases. This coupling is controlled by the pressure equilibrium, characterised by very small time scales. The relaxation parameter μ quantifies the mechanical equilibrium time scale. The source terms in the specific internal energy equations appear to represent the exchange of the specific internal energies due to the pressure work. The mixture speed of sound in this model is computed from: where Y denotes the mass fraction. As Eq. 1 denotes, the effects of viscosity, heat conductivity, surface tension, and phase transition are neglected in the present study. which presents a monotonic variation with volume fraction [45]. The numerical solution of the Eq. 1 is obtained through three major steps at each temporal loop: 1. Solving the hyperbolic part of the system. In this step, the pressure disequilibrium is assumed and the relaxation terms are ignored. This treatment gives a hyperbolic system for the conservative variables that is solved using an approximate Riemann solver with a finite volume scheme. 2. Converging to an equilibrium pressure by solving the relaxation system described in B. 3. Correcting of the solution by enforcing the total energy conservation.

Hyperbolic step
Considering the homogenous part of the PDEs in this step, Eq. 1 is solved using a finite volume Godunov method [58] with the secondorder MUSCL scheme [59] employed to reconstruct the primitive variables at the cell boundary. Moreover, the HLLC approximate solver [60] is adopted to solve the Riemann problem at each cell boundary as an appropriate choice for the present method [45,61,55].
Using the computed inter-cell fluxes, the solution of the conservative and non-conservative variables can be evolved on the entire time step. The conventional Godunov scheme to update the conservative part of the system reads: in which: The subscripts (i,j) stands for the finite volume cell index in (r,z) direction and superscript n shows the time step. Superscript ' * ' denotes the perturbated state. Calculation of F * cons and G * cons using the HLLC Riemann solver is explained in [45] and provided in C of the present work. The non-conservative part of the equations is updated following by approximating the volume integral with a midpoint rule and the divergences with a centred scheme [45]. Using this approximation, the volume fraction is updated as: Assuming that the product (αp) n i,j is constant during the time step, the non-conservative internal energy equations can be approximated as [45]: The approximation of the internal energy at this step is not crucial as it is used only to estimate the phasic pressures which will be corrected later in the relaxation step [45]. Finally, the solution is advanced using a two-step time integration [32]: where L R contains the right hand side of Eqs. (4)-(6).

Relaxation step
The hyperbolic step leads to different phasic pressures whereas the numerical solution should converge to a unique pressure to fulfil the mechanical equilibrium in the interface. This is achieved through the relaxation step. The concept behind the relaxation approach stems from two observations: firstly, when a rarefaction or shock wave passes through two phases having different pressures, the volume of each phase must change in order that pressures tend to equilibrium. The first relaxation term in s rlx of Eq. 1 represents this volume fraction expansion with rate μ. Secondly, a pressure work is associated with the volume change of the phases. This is reflected by the last two relaxation terms in s rlx of Eq. 1. Physically, μ depends on the mechanical properties of the fluids as well as the mixture topology [62,63]. The stiff pressure relaxation with μ→∞ results in instantaneous pressure equilibrium at the interface at any time [62,36]. In [45,56], it is demonstrated the instantaneous pressure equilibrium is a valid assumption. From the relaxation system described in B, the following equations can be derived: shows the specific volume of phase k. The interfacial pressures p I in both phases are considered to be equal such that the mixture energy is conserved. A possible estimation of p I is the mixture pressure computed after the relaxation step [45,56], which fulfils the entropy inequality (more details can be found in [64]). This allows the construction of a non-linear algebraic system: where superscripts (b) and (a) denote the values before and after the relaxation step, respectively. Considering the five unknowns e (a) , three more equations are needed to close the system. First, since (αρ) k is conserved during the relaxation, the saturation constraint ∑ k α k = 1 reads: The two more required equations are extracted from the phasic equations of state, which express the internal energy of the phase based on its pressure and density. This is straightforward when using simplistic equations of state, such as the ideal gas and stiffened gas EoSs described in A in their parametric forms, which read as: In this specific case, an analytical solution for the system of Eqs. (11)-(14) exists, which is described in [45,65]. However, in the case of complex equations of state, there is no analytical solution for this system. Therefore, an iterative solution should be tailored based on the particular formula of the utilised EoS. Herein, we present a general algorithm for this system based on tabulated data. In this regard, residuals associated with Eqs. (11)-(13) are defined: This system can be solved iteratively using the multivariable Newton's method. As the pressure and temperature are the inputs in the tabulated data, the system is solving for T (a) 1 , T (a) 2 , p (a) through the following steps:    4. The following iterative system in constructed based on the Newton's method: where J − 1 stands for the inverse of the Jacobian matrix:

5.
To increase the convergence speed and avoid instabilities, an underrelaxation coefficient c urf = 0.5 was applied at each loop as: where superscript (a) (n+1) ☆ denotes the updated values at the end of each iteration of the Newton's method. The following criteria based on the change of the variables over the iterations is considered at each point: where ε = 10 − 2 is sufficient for convergence to be achieved. 6. Calculating the equilibrium pressure and the corresponding phasic densities, the volume fractions can be updated as (αρ) k is conserved for phase k. From the relaxation step, the phasic densities and volume fractions are estimated properly as inferred from the numerical tests of [66].

Re-initialisation
After the hyperbolic step, the non-conservative internal energies do not infer a common relaxation pressure. In addition, the phasic internal energies, solved in a non-conservative manner, do not necessarily correspond to the total energy. Thus, a re-initialisation step is introduced to ensure that the non-conservative internal energies are consistent with a common equilibrium pressure and correspond to the total energy. This is summarized in the following steps: 1. The mixture internal energy e m is extracted from the total and kinetic energy computed in the hyperbolic step: 2. The mixture rule for the internal energies is considered as: in which the left hand side e m is known from the previous step while Y 1 and Y 2 are available from the relaxation step. Therefore, with the substitution of e 1 and e 2 with the corresponding pressure and densities based on the equations of state, the pressure will be the only unknown. In the case of the stiffened gas EoS, Eq. A.1, we obtain the following equation: For more complex equations of state, however, there is no analytical solution of Eq. 21. Similar to generalisation made in the relaxation system, we present an iterative method valid for any arbitrary equation of state. In this regard, the residual based on Eq. 20 is considered: Newton's method for pressure reads: where ∊ ′ p is the partial derivative of the residual function with respect to the pressure estimated as: where Δp represents a small change in pressure and can be estimated based on the pressure from the previous loop Δp = ξ p p for which ξ p = 10 − 3 is recommended. Moreover, the initial guess values are considered based on the relaxation step. Similar to the relaxation step, an under-relaxation treatment is also considered to ensure stability. The pressure is computed when the solution converges ⃒ <ε; a suggested value of ε = 10 − 3 has been used.
3. Finally, the internal energies are recomputed with the pressure obtained in step 2. It is now ensured that the updated internal energies are in agreement with the total energy conservation.

Results and discussion
In this part, three spherical collapse cases with different initial pressure ratios are first presented, followed by the non-spherical collapse of a bubble near a rigid wall excited by a pressure pulse corresponding to conditions similar to those generated by commercial lithotripter ultrasound systems. Moreover, the effect of the distance between the bubble and the rigid wall is considered. A cavitation case is also presented in D as a benchmark test. In all simulation, the gas phase is assumed to be non-condensable and the effects of viscosity, heat conductivity, surface tension, and phase transition are neglected.
Each phase contains a small volume fraction of the opposite phase α min = 10 − 6 in the initial setups to ensure the hyperbolicity of the system, as recommended in [45]. Moreover, the monotonized central slop limiter is used for the MUSCL reconstruction scheme explained in [45]. The time step is varying based on the CFL number which is set to 0.5. The first case considered here is the symmetric collapse of an isolated air bubble with initial radius R 0 = 1 mm submerged in infinite water at rest; the different investigated are indicated in Table 1. This test is performed in 1D spherical coordinates. Initially, the pressure inside the air bubble p air is uniform in r = [0, R 0 ] while the surrounding pressure in r = (R 0 , L D ] increases gradually towards the far-field p f [1]: where L D indicates the domain size. The grid cells are uniformly distributed in two regions with different resolutions. In the first region r = [0, 3R 0 ] containing the bubble and its neighbourhood, 3N R0 cells are uniformly placed where N R0 denotes the number of cells per initial radius R 0 . The number of N R0 = 100 is sufficient to obtain the converged   = 2N R0 , 6N R0 , 10N R0 cells for case 1, 2, and 3, respectively, to be large enough to avoid any possible interaction between the wave reflection from the far-field and the bubble. Reflective and transmissive boundary conditions [60] are used for the bubble centre and the far-field region, respectively. Results obtained with the ideal and real gas equations of state are also compared to those obtained with the Keller-Miksis model. It should be mentioned that the Helmholtz EoS failed in the simulation of case 1 as the bubble properties in this case exceed the valid range of the Helmholtz EoS described in Section 2. To make the comparisons more clear, we use initial radius and Rayleigh collapse time to nondimensionalise the radius and time respectively as follows [33]: As can be seen in Figs. 2a, 2c and 2e depicting cases 1, 2 and 3 respectively, the bubble undergoes a compression, the rate of which depends on the initial pressure ratio, followed by a rebound. The comparison with the Keller-Miksis model shows that the present method captures the compression and expansion rate with satisfactory accuracy for all cases. The temporal change of the space-averaged gas temperature inside the bubble is plotted for all cases in Figs. 2b, 2d, 2f using both the real and ideal gas EoSs. It is found that the predicted temperatures obtained with the three real gas EoSs are very similar, as expected considering the results presented previously in Fig. 1. On the other hand, the difference between the temperatures predicted by the real gas EoSs and ideal gas EoS is significantly affected by the initial pressure ratio. It is observed that for the violent collapse, case 1, the space-averaged temperatures obtained with the real-gas EoSs are ≈ 33% lower than the value predicted by the ideal gas EoS. The difference is negligible, however, in case 3 where the collapse is mild. The maximum temperature achieved during bubble collapse is reported for all cases investigated in Table 2.  To get more insight into the bubble spatial temperature variation, the spatio-temporal distribution of the gas temperature in r = [0, R 0 ] obtained with the ideal gas and RKPR EoSs for case 1 is depicted in Fig. 3. In this figure, the vertical axis shows the non-dimensional space inside the bubble while the horizontal axis denotes the non-dimensional time. This representation in useful to illustrate the temperature locally inside the bubble during the entire simulation. It is evident that the bubble collapse undergoes a nearly isothermal process in the initial stage due to the slow collapse rate, followed by adiabatic heating. As the bubble approaches to its minimum size during the collapse, the temperature rises due to the very high compression rate and reaches a local maximum values of ≈ 10, 000 K and ≈ 6, 000 K predicted by the ideal gas and the RKPR EoS respectively. Subsequently, the bubble cools down during the expansion phase where gradients of temperature are observed.
It should be noted that the number of iterations required for the relaxation to converge varies in time. To demonstrate this, the average of this value in the whole domain has been reported in Fig. 4 for case 1 with the RKPR EoS. Accordingly, the required number of iterations reaches its maximum during the collapse to reach the convergence. The CPU time of the serial computation for this simulation is 320.85 s on an Intel Core i7-8850U CPU @1.8 GHz.

Shock-induced non-spherical collapse close to a rigid wall
Moving towards a case with more practical interest, a shock-induced bubble collapse near a rigid wall surface is examined. The bubble-wall arrangement resembles that of a lithotripter system. This test case was first introduced in the work of [50], where they studied the wall pressure subjected to the bubble collapse. In this setup, infinite impedance for the kidney stone is assumed to avoid any wave absorption in the boundary. It was shown that the wall pressure reaches values of GPa depending on the initial stand-off distance as well as the pulse width and amplitude. Herein, the focus is on the collapse temperature with the ideal and RKPR EoSs. Moreover, the wall pressure is depicted for various initial stand-off distances.
The compressive shock front from the upper boundary demonstrated in Fig. 5a represents the lithotripter pulse without the tensile part propagating in time; this is based on an analytical function described in [50]. Initially, the pressure is atmospheric in the whole domain and the water and air densities are ρ water = 998.2 kg/m 3 and ρ air = 1.125 kg/m 3 , respectively. To reduce the computational cost, the case is simulated in 2D axisymmetric coordinates instead of the full 3D configuration. A   Fig. 10. Gas temperature variation over time in the 2D axisymmetric simulation of the non-spherical shock-induced collapse with the real gas (left column) and ideal gas (right column) EoS for H 0 = 2R 0 . a) t * ≈ 11.54, b) t * ≈ 12.82, c) t * ≈ 13.65, d) t * ≈ 14.11. schematic of the geometry and the mesh is presented in Fig. 5b. Different resolutions in r = [0, 1.2R 0 ], r = (1.2R 0 , 2R 0 ] , and r = (2R 0 , 4R 0 ] are used in the r-direction with the total number of N r = 400 cells. In the zdirection, however, the grid is uniform in z = [0, 1.2R 0 ] while the Vinokur function [67] is used in z = (1.2R 0 , 15R 0 ] for grid stretching with the total number of N z = 750 cells. The grid independence study is provided in E, showing that considering N R0 = 150 cells in the initial setup is sufficient to obtain grid independent solutions. The reflective boundary condition is used for the axis of symmetry whereas for the right side and the bottom wall, the non-reflective and no-slip boundary conditions have been used, respectively. The bubble has an initial radius of R 0 = 0.05 mm while the initial stand-off distance case is 2R 0 .
To be consistent with the reference [50], the variation of the bubble volume normalised with its initial value V * = V/V 0 over nondimensional time t * = tc L /R 0 is plotted in Fig. 6 where c L = 1647 m/s is the reference speed of sound. The results obtained with the ideal and real gas EoSs are compared with the study of [50]; overall, good agreement is achieved. It can be seen that the dynamics of the bubble is not affected by the gas thermodynamics. On the other hand, it was seen in the spherical bubble collapse cases that the bubble temperature predicted by the three real gas EoSs was nearly the same as shown in Fig. 2. Therefore, the RKPR EoS is used as the real gas EoS for the rest of the results.
In Fig. 7a, it is observed that the bubble compression results to bubble space-averaged temperatures up to 6, 500 K when the ideal gas EoS is utilised. The temperature predicted with the RKPR EoS is approximately 3, 800 K which is ≈ 41% lower than the ideal gas EoS prediction. This noticeable temperature difference at p b ≈ 4 GPa, which is corresponding to the compression ratio of 4 × 10 4 , is in agreement with the comparison made in Section 2, Fig. 1 where ≈ 35% temperature difference between the ideal and real gas prediction for an isentropic compression at similar compression ratio was observed. Fig. 7b, however, shows a less significant difference in the bubble pressure which is only ≈ 2.3% based on our data. Similarly, in Fig. 7c, it is observed that the predicted wall pressure averaged in r = [0, R 0 ] is less affected by the gas EoS due to the sufficiently large stand-off distance.
Simulations performed with the RKPR EoS for different stand-off distances namely H 0 = 1.1R 0 , H 0 = 2R 0 , and H 0 = 3R 0 as shown in Fig. 8a-8c. It is observed that when the stand-off distance is minimum, i. e., H 0 = 1.1R 0 , the collapse forms in a more asymmetric shape compared to the two other cases. Therefore, a less amount of energy can be concentrated inside the bubble. As a result, the maximum bubble temperature and pressure are lower when H 0 = 1.1R 0 . On the other hand, the wall pressure peak, averaged in z = [0, R 0 ], is the highest in this case as the shock immediately hits the wall.
Numerical Schlieren [68] is used as a useful gradient-based function for visualisation of the formed waves as well as the interface location [35,69]: where ψ is a scaling parameter to improve the visibility of the waves for which a value of 50 is used in the simulations. In the contour plots of Fig. 9, the temporal evolution of the pressure field p and velocity vectors (on the left-half), and the numerical Schlieren (on the right-half) are illustrated at the selected times, namely a) t * ≈ 11.54, b) t * ≈ 12.82, c) t * ≈ 13.65, d) t * ≈ 14.11 corresponding to different collapse stages. Moreover, the line plots of Fig. 9 represent the wall pressure p w at each time. The simulations are obtained for H 0 = 2R 0 with the RKPR EoS while no substantial difference for this set of variables were observed in the case of the ideal gas EoS. It can be seen that the emitted pressure pulse from the lithotripter hits the bubble at the top boundary creating a reflected rarefaction wave and a transmitted shock wave. The induced pressure gradient onsets the bubble collapse, as seen in Fig. 9a (t * ≈ 11.54). The transmitted wave is then reflected from the rigid wall and hits back the bubble from the lower boundary, compressing it further. At this stage, the peak pressure of the wall is ≈ 0.055 GPa. The liquid jet formation is evident in Fig. 9b (t * ≈ 12.82) as well as an emitted shock wave due to bubble collapse. Thirdly, the bubble expands into a toroidal-like shape and the emitted shock wave from the collapse travels toward the rigid wall, as seen in Fig. 9c (t * ≈ 13.65). Up to this point, the wall pressure peak is still less than 0.08 GPa. Lastly, the bubble expands further and the shock wave hits the rigid wall and abruptly increases its pressure to the peak value of above 0.45 GPa, as depicted in Fig. 9d (t * ≈ 14.11). This shows that the peak value of the wall pressure is strongly influenced by the collapse shock wave. Also, observing the wall pressure above 0.37 GPa in |r| < 30 μm it seems that the wall can retain the high pressure value as the shock is passing outward from the centre. While the spatial average of the gas temperature presented in 7 is appropriate to monitor the bubble temperature over time, we are able to extract more insight about the temperature distribution and the real gas effects locally at different collapse stages with a 2D representation. Therefore, the evolution of the gas temperature is shown in Fig. 10 at the same selected times considered for Fig. 9, i.e., t * ≈ 11.54, b) t * ≈ 12.82, c) t * ≈ 13.65, d) t * ≈ 14.11. The results obtained with the RKPR EoS (left column) are compared with those with the ideal gas EoS (right column). In both cases, the temperature distribution inside the bubble is inhomogeneous as expected in non-spherical collapse. This temperature first increases due to the adiabatic during the compression phase. At this stage, the difference between the predictions by the ideal and RKPR EoSs is minor, with the maximum of ≈ 60 K in the spot shown in Figs. 10a (t * ≈ 11.54). However, this difference becomes substantial at the moment of collapse. At the selected t * ≈ 12.82, which is very close to the collapse moment, the maximum temperature predicted by the ideal gas EoS reaches ≈ 10000 K in the very small red spot in the centre as shown in Fig. 10b R which makes a difference of ≈ 3800 K with the maximum prediction by the RKPR EoS in Fig. 10bL. Although this maximum temperature difference happens only in a very small spot, the RKPR EoS predicts lower collapse temperature even in the surrounding of the spot. As the bubble expands, cooling of the gaseous content is expected. This is clearly seen in Figs. 10c (t * ≈ 13.65) and 10d (t * ≈ 14.11) with the maximum temperature difference of ≈ 300 K and ≈ 200 K, respectively.

Conclusion
A numerical model extending the work of [45] to account real-fluid EoS via tabulated data has been developed and applied to gas collapse cases. Both spherical and non-spherical collapse cases have been considered and compared with simulations obtained with the ideal gas EoS. The real gas effects become more dominant when the collapse is more violent, leading to a ≈ 33% difference (corresponding to ≈ 3050 K) in the space-averaged spherical collapse temperature. Also, a difference of ≈ 41% (corresponding to ≈ 2700 K) is observed for the spaceaveraged non-spherical collapse temperature. It is also shown that the difference is even more evident if the local extreme collapse temperature is considered. Therefore, it is concluded that the ideal gas assumption fails for temperature prediction regardless of the collapse sphericity.

CRediT authorship contribution statement
in which subscript 0 denotes the ideal state.
The following parameters are used in Table A.3: The specific heat capacities at constant volume and pressure are computed: Moreover, the speed of sound can be derived as: where s stands for the specific entropy.

A.2.3. c) Helmholtz energy EoS
The Helmholtz energy EoS is based on an ideal and residual Helmholtz energy formulation, where the residual term α r is a large polynomial/ exponential function, calibrated with experimental data [76]:
Moreover, δ = ρ ρ j is the reduced density and τ = Tj T is the reciprocal reduced temperature where subscript j denotes the value at the maxcondentherm. The residual Helmholtz energy contribution can be computed as: for which the coefficients N k and values of i k , j k , and j k can be found in [76] as well as the ideal gas contribution.

Appendix B. Relaxation system
The relaxation step contains the following equations [45]:

Appendix C. HLLC Riemann solver
The HLLC solver is being used extensively calculating the fluxes with high accuracy. This solver is recommended for the present method [45,55,65]. In this solver, the bounds for the minimum and maximum signal velocities present in the solution of the Riemann problem can be estimated as: where c stands for the mixture sound speed. In case of using the parametric SG equation of state, the sound speed for each phase k reads: When the tabulated equation of state is considered, the speed of sound is computed based on p k and ρ k through an interpolation using the tabulated data. The intermediate wave speed can be calculated using the HLL approximation using the mixture density and pressure: In the absence of the relaxation effects, the volume fraction is constant along the fluid trajectories: α * kR = α kR , α * kL = α kL . (C.12) Since the volume fraction does not change across the left and right wave, the density reads: To compute the internal energy jumps the Hugoniot relation is required: If the parametric equations of state are considered, the internal energy can be expressed as a function of pressure and density: .
(C.15) Therefore, there will be only one unknown in the Eq. (C.14). In case of using the SG equation of state, for instance, using Eq. (A.1): Subsequently, the internal energy e * k can be obtained as a function of the pressure and density using the equation of state. However, the Hugoniot relation in Eq. (C.14) should be solved iteratively for energy and pressure in the case of using the tabulated EoS for the gas since there is no analytical expression relating the pressure and the internal energy. Therefore, the residual is defined based on Eq. (C.14) and the Newton method is used to minimise the error iteratively.

D.1. Cavitation
Cavitating of a water-air mixture is simulated with the RKPR EoS and IG EoS and the results are compared with the exact solution. Initially, the pressure is p = 1 bar everywhere in the domain x = [0, 1] m while a velocity discontinuity is placed in the middle of the domain x = 0.5 m. The flow velocity on the right side has u = 100 m/s while it is set to u = − 100 m/s on the left side. The tube is filled with the liquid water with ρ = 1000 kg/m 3 while there is a very small air volume fraction α a = 10 − 2 in the domain.
Existence of the gas with the left and right expansion waves from the centre generate vacuum in the middle of the domain where the liquid pressure drops, and a cavitation pocket is formed. Good agreement between the results with ideal gas EoS and the exact solution at t = 1.85 ms from [45] is demonstrated in Fig. 11 in which the two interfaces are captured with a high accuracy with 1000 cells. It should be noted that real gas effects are minor in this regime of pressure. Therefore, it is observed that the results with RKPR EoS is similar to the results with the ideal gas EoS.

Appendix E. Grid convergence
The discrepancy between the numerical results and the Keller-Miksis solution for the low initial pressure ratio case can be served as the error when the grid resolution is analysed [33]. In this regard, the error is defined as the L 2 error for case 3: in which ε is the error computed in the time [0, 2t c ] in N t time steps where t c is the collapse time and R(t i ) and R KM (t i ) are the radius from the simulation and the Keller-Miksis model, respectively. The results shown in Fig. 12a represent a first order convergence which is the expected rate in the case of the discontinuous flows [60]. The grid resolution analysis for the same case is shown in Fig. 12b. Grid convergence for the shock-induced non-spherical bubble collapse is also shown in Fig. 13.