A smoothed particle hydrodynamics study of the collapse for a cylindrical cavity

In this study, we propose a mesh-free (particle-based) Smoothed Particle Hydrodynamics model for simulating a Rayleigh collapse. Both empty and gas cavities are investigates and the role of heat diffusion is also accounted for. The system behaves very differently according to the ratio between the characteristic time of collapse and the characteristic time of thermal diffusion. This study identifies five different possible behaviours that range from isothermal to adiabatic.


Introduction
The term "cavitation" describes a phenomenon composed by two distinct phases: first, a vapour cavity, also called vapour bubble or void, develops and rapidly grows in a liquid phase; subsequently, the vapour cavity rapidly collapses generating strong shock waves.
Cavitation causes erosion and it is mostly undesirable in engineering applications such as turbo-machines, propellers, and fuel injectors [1][2][3]. However, other applications such as ultrasonic cleaning or cataract surgery [4][5][6] are specifically designed to take advantage of the erosion power of the collapsing bubble.
According to the circumstances, the bubble collapse can follow two distinct, but similar, mechanisms [7] called, respectively, Rayleigh collapse and shock-induced collapse. During the Rayleigh collapse, the collapse is driven by the pressure difference between the surrounding liquid and the cavity. In this case, if the pressure field is perfectly isotropic, the bubble maintains a spherical shape during the whole duration of the collapse. Shock-induced collapse is caused by the passages of a shock-wave through the bubble. In this case, the spherical shape is not preserved and the bubble folds in the shock direction.
Our current understanding of cavitation is based on three different approaches: (i) theoretical investigations, (ii) experiments and (iii) computer simulations.
The first analytical study of an empty cavity surrounded by an incompressible fluid at given pressure was carried out by W. H. Besant (1859) [8], who obtained an integral expression for determining the time required for the cavity to collapse due to the effect of a constant external pressure. Sixty years later, Lord Rayleigh [9] was able to integrate this equation determining that, during the collapse, the pressure of the liquid near the boundary exceeds the pressure of surrounding liquid. Plesset [10] introduced the effect of surface tension and viscosity obtaining nature [37,38]. The method has been validated for wide range of applications such as explosion [39], underwater explosion [40], shock waves [28, 41,42], high (or hyper) velocity impact [43], water/soil-suspension flows [44], free surface flows [45,46], nano-fluid flows [47], thermo-fluid application [48]. Moreover, SPH is also a component of the Discrete multi-physics simulations [49][50][51][52][53] SPH bases its discrete approximation of a continuum medium on the expression where f(r) is any continuum function depending of the three-dimensional position vector r, while W is the smoothing function or kernel. The kernel function W defines the extension of the support domain, the consistency, and accuracy of the particle approximation [25]. When the computational domain is divided in computational particles with their own mass, m = ρdr, it is possible to rewrite Eq 1 in particle form where m i , ρ i and r i are mass, density and position of the i th particle. Within the SPH framework, it is possible to use Eq 2 to discretise a set of equations such as the continuity equation which in particle form becomes the momentum equation which in particle form becomes where P ij is called artificial viscosity and was introduced by Monaghan [41] for simulating shock waves, having the expression where β is the dimensionless dissipation factor, c i and c j are the speed of sound of particle i and j.
The energy conservation equation which in particle form becomes where κ is the thermal conductivity. The first term of the right side of Eq 9 is the particle form of the sum of the reversible rate of the internal energy increase by compression and the irreversible rate of internal energy increase by viscous dissipation; the second term is the particle form of the rate of internal energy increment by heat conduction following Fourier's Law. The thermal conductivity is related to the thermal diffusivity, α, by the relationship where c p is the specific heat capacity.

Kernels function
In this work, two kernels (Lucy Kernel function and quintic spline) are used and their effect on the accuracy of the results compared. The Lucy kernel function [33] Wðq; hÞ ¼ is one of the simplest kernels used in literature. Where q = |r − r 0 |/h, s is a parameter used to normalise the kernel function, which, for one, two and three dimensional space is, respectively, 4h 5 , ph 2 5 and 16ph 3 105 . The quintic spline is a piecewise kernel function [54] WðQ; hÞ ¼ s ð3 À qÞ 5 À 6ð1 À qÞ 5 þ 15ð1 À qÞ 5 ; 0 < q < 1 ð3 À qÞ 5 À 6ð1 À qÞ 5 ; where q = |r − r 0 |/h and s for one, two and three dimensional space is, respectively 1 120h , 7 478ph 2 and 3 359ph 3 . In the quintic kernel is normally more accurate, but at the expenses of higher computational costs because it requires a neighbor list three times larger than the Lucy kernel [25].

Model
Problem description. In the Rayleigh collapse, the driver force is the pressure difference between the pressure in the liquid, p 1 = P L , and the pressure in the cavity, p b .
Two different scenarios are analysed: empty cavity collapse, where the cavity is void, with p b = 0, surrounded by a liquid phase, and vapour cavity collapse, where the cavity is filled of a non condensable gas with an initial pressure equal to the vapour pressure of water at the temperature T 0 , p b = p sat (T 0 ) and density equal to the density of an ideal gas at that pressure and temperature, ρ b (p b , T 0 ).
The liquid phase is water with ρ L = 1000 kg/m 3 , P L = 5 MPa and T 0 = 300 K. The liquid pressure of 5 MPa has been chosen as it is commonly reached in various hydraulic applications [55]. Given the liquid temperature, the pressure in the bubble is p b = 3.55 kPa with a density of ρ b (p b , T 0 ) = 2.7 � 10 −3 kg/m 3 . At the short timescale considered, water is compressible. The equation of state for compressible water is discussed later on.
In rapid collapse, the water vapour is considered trapped within the liquid, assuming zero mass transport across the interface. This simplification is justified because mass transport mechanisms across the interface and non-equilibrium condensation require higher timescales to play an effective role in the collapse phase [56]. The short timescale also justifies neglecting the surface tension in modelling the collapse.
The short timescale of the phenomenon may suggest an adiabatic collapse [20,57]. However, especially in the last stage of the collapse, the high temperature gradient between gas and liquid phases could introduce a non-negligible heat transfer between the two phases [21,58]. In this study, both scenarios (e.g. adiabatic and non-adiabatic collapse) are investigated.
SPH model. The axisymmetric water domain is shown in Fig 1. The domain is dived in three concentric regions, delimited by three different radiuses, where different types of computational particles are used.
Cavity (r < R 0 ): inside this region particles are removed (in the case of empty cavity) or modelled as non-condensable gas following a gas phase equation of state (EOS) in the case of vapour collapse. In the rest of the paper, particles inside the cavity (when present) will be referred as particle Type 1.
Liquid (R 0 < r < R S ): inside this region particles are modelled as compressible fluid following a liquid EOS. Particles inside the liquid will be referred as particle Type 2.
Shell (r > R S ): inside this region particles are modelled as fluid with a fixed position and density to represent the boundary conditions of the system and maintain a fixed pressure at the boundaries. Particles inside the shell will be referred as particle Type 3. We also run several simulations with cubic control volumes and periodic conditions that account only for Type 1 and Type 2 particles. The results do not change and, therefore, we prefer the system in Fig 1 that overall requires less computational particles.
To avoid compenetration between gas and liquid particles during the gas cavity collapse, we employed a penalty force, similar to the one used by Liu et al. [40] between these types of particles.
with C = 10 −4 , γ = 9 and σ equal to the initial particle spacing. The initial radius of the cavity is R 0 = 100μm (typical radius of a collapsing cavity [26,59]). The ratio R C /R 0 = 30 is used as a compromise between computational cost and accuracy. Different R S has been tested, as explained in the Hydrodynamic section.

Equation of state.
To solve the set of Eqs 3-8, an EOS that links pressure P and density ρ is required. Each phase requires a different EOS: in this work multiple EOS are used and compared.
Liquid EOS. For liquids, we used and compared two EOS: the Tait and the Mie-Gruneisen EOS.
The Tait equation is probably the most used EOS in SPH to model water where c 0 is speed of sound of the liquid and ρ 0 is the reference density. The Tait EOS takes in account the compressibility of the liquid. Is possible to regulate the compressibility by selecting the appropriate sound of speed [60] in Eq 14.
For simulating underwater explosion Liu et all [40] used the Mie-Gruneisen EOS [61] to model the water as a compressible fluid having different expressions for compression and expansion state. Shin et al. [62] derived a polynomial expression for both compression and expansion states: for compression state, while for expansion state, Where μ = ρ/ρ 0 − 1 and e is the specific internal energy. The coefficients are a 1 = 2.19 � 10 9 N/m 2 , a 2 = 9.224 � 10 9 N/m 2 ,a 3 = 8.767 � 10 9 N/m 2 , b 0 = 0.4934 and b 1 = 1.3937 evaluated for water with ρ 0 = 1000 kg/m 3 and c 0 = 1480 m/s. Vapour EOS. The vapour phase in the cavity is modelled as a non-condensable gas. In our simulations, we used and compared two EOS: the ideal gas EOS and the NASG EOS.
The ideal gas EOS, for the pressure, is given by

PLOS ONE
temperature where γ = c p /c v is the capacity heat ratio, M m the molar mass of the gas and R is the ideal gas constant. The NASG EOS, which is a multiphase EOS is discussed in the next section. Multiphase EOS. Le Métayer & Saurel [63] combined the "Noble-Abel" and the "Stiffened-Gas" EOS proposing a EOS called Noble-Abel Stiffened-Gas (NASG), suitable for mulfiphase flow. The expression of the EOS does not change with the phase considered, and, for each phases, is possible to determine both the pressure and temperature as function of density and specific internal energy. Pressure-wise the expression of NASG is and temperature wise where P, ρ, e, and q are, respectively, the pressure, the density, the specific internal energy, and the heat bond of the corresponding phase. γ, P 1 , q, and b are constant coefficients that defines the thermodynamic properties of the fluid. The coefficients for liquid water and steam used in our simulations are given in Table 1. Software for simulation, visualisation and post-process. All the simulation were run with the open source code simulator LAMMPS [64,65]. Visualisation and data post-processing were generated with the Open Source code OVITO [66].

Hydrodynamic
Empty cavity. Different simulations have been run to assess the quality of results with respect of numerical parameters such as number of computational particles, kernel function and time step. Preliminary simulations have been run using both Lucy (Eq 11) and quintic spline (Eq 12) kernel functions obtaining similar results. Therefore, we chose the Lucy kernel over the over the quintic because it requires less computational cost because accounts for a smaller neighbour list. In all cases smoothing length and the dissipation factor are h = 1.3 � dL, where dL is the initial particle spacing, and β = 1, coherent with literature in shock-wave problems [25,42].  Note that our dimensionless radius does not goes to zero, but it rebounds, like the axisymmetric Raylerigh-Plesset equation, this is explainable with the particle nature of the SPH method: with particle methods, in fact, there is always a small spacing between particles.

PLOS ONE
Based on the analysis of this section we decided, for an empty cavity collapse, to use the parameters shown in Table 2.
All the parameter listed are chosen as the best compromise between speed an accuracy.  Where e 0 is the initial internal energy of particles. The timestep t s = 10 −12 s was chosen since it is the largest timestep where the internal energy decreases after the collapse as required by the physics of the problem. The parameters used are summarised in Table 3.
Comparison with the axisymmetric Rayleigh-Plesset equation. The Hydrodynamic of the model is compared with the solution of the axisymmetric Rayleigh-Plesset (ARP) Eq [7] for both the empty and vapour cavity. Fig 4 shows the dimensionless radius, R(t)/R 0 , plotted against dimensionless time, t/t c , of our model against the solution of equation ARP for the empty collapse. R 0 is the initial radius of the cavity, t c = 2.76μs is the collapsing time obtained with the ARP.
In our model, the cavity collapse slightly faster than the theoretical, leading to t/t c � 0.98 instead of 1. This difference is explainable with the compressibility of the liquid [13]. The theoretical model assumes that water is perfectly incompressible, while the Tait EOS in the SPH model accounts for the compressibility of water. At these timescales, the compressibility cannot be neglected, and, from this point of view, our compressible SPH model should be more accurate than the theoretical, fully incompressible model. It is also important to highlight that, in our simulation, the parameter c 0 (sound speed in the medium) in the Tait EOS, we use is 1484 m/s, which is the actual speed of sound in water at 25˚C.
According to our calculations, the effect of the compressibility affects the rate of collapse. At the beginning, the compressibility produces a higher collapsing rate because a compressible fluid fills the void in the cavity faster than an incompressible fluid. As the cavity shrinks, however, the curvature of the cavity acts as an arch and the speed of the collapse slows down more than in the more rigid (incompressible) case. Overall, these two effects cancels each other out   The considerations done for the empty collapse are still valid for the vapour cavity case. However, additional discussion is required for the final phase of the collapse and the rebound phases: unlike the empty cavity (see Fig 6a) when the vapour cavity approaches the final phase of the collapse, the cavity loses the cylindrical symmetry (Fig 6b), and differs from the ARP.

PLOS ONE
This "artificial" asymmetry in the rebound phase is attributable to low resolution occurring when, during the last stage of the collapse, the size of the caivty is comparable to the size of the smoothing length. However, this work focuses only on the bubble collapse (as usual in computer simulations of cavitation e.g. [7,26,59]) and the rebound phase is not considered.

Pressure field
Pressure field in the liquid (empty cavity). Initially, at t = 0, the pressure is uniform along the domain. As the cavity shrinks, due to the pressure difference between the liquid and the cavity, the liquid starts to fill the cavity. This causes a decrement in pressure in the liquid generating a low-pressure wave that moves through the liquid phase (see Fig 7).
As the collapse proceed, a high-pressure area arises near the cavity border (see Fig 8a and  8b) that abruptly increases reaching the max at the collapse (see Fig 8c and 8d). Locally, the max pressure calculated is around 120 MPa. This value is one order of magnitude lower than the theoretical value calculated by Hickling and Plesset [68] for the 3D collapse, but this difference is consistent with the fact that our model refers to a 2D collapse [67,69]. This is also reflected in the difference in the collapse time between 2D and 3D [9,70]

PLOS ONE
The collapse generates a high-pressure wave (Fig 9 and 9b) that moves away form the cavity (Fig 9c and 9d). There is a theoretical reason for the hexagonal patterns in Figs 8 and 9, which is discussed in the next section.
With the absence of heat diffusivity the presence of the vapour in the cavity does not affect significantly the pressure in the liquid and, therefore, pressure fields for the vapour cavity collapse are not shown here.
Acoustic diffraction. As mentioned in the previous section, when the empty cavity reaches the minimum radius, a high-pressure shock wave is generated and propagates in the liquid phase. After the wave bounces back, it loses its spherical symmetry and assumes an unphysical hexagonal symmetry. This is a numerical artefact and depends on the fact that, below a certain raito R/dL, the initial particle resolution is not adequate to correctly describe the cavity shape (the reader can compare Fig 7, where R/dL = 19.27, with Fig 8, where R/dL = 4.50-1.75). The cavity assumes a hexagonal shape (see Fig 8a), which is related with initial hexagonal particle distribution of the model. When the high-pressure shock wave bounces back, therefore, it propagates from a hexagonal cavity rather than a circular one.
This behaviour closely resembles light diffraction from a hexagonal aperture, Fig 10. Light diffraction, in fact, follows specific patterns [71] defined by the shape of the aperture. Fig 10 shows the similarity between the light intensity pattern generated by diffraction trough a hexagonal opening and the pressure intensity pattern of the wave generated at the

PLOS ONE
collapse of the cavity. The pressure peaks and valleys in Fig 10b (and Fig 9a), therefore, are the results of Fresnel like positive and negative interference of the interfering diffracted waves rather than the result of effect of numerical instability.
This issue, however, only occurs at the end of the collapse, when the size of the cavity is comparable to the smoothing length and does not affect the overall collapsing time.

Thermal dffects
During the collapse, the compression of gas in the bubble generates heat. This heat, in turn, can affect the dynamic of the collapsing bubble [14]. When thermal effects are absent, or negligible, the collapse is "inertially controlled" as in the previous Section. When the thermal effects are not negligible, the collapse can be "thermally controlled". In a thermally controlled collapse, the bubble dynamic differs from the inertially controlled because the thermal terms in the ARP equation are not negligible.
Two scenarios are analysed. In the adiabatic collapse, only the first term of Eq 9 is accounted for. In the heat diffusive collapse, both term are enabled to model heat transfer between gas and liquid.
Finally, the temperature peak in the gas cavity is investigated in relation to the ratio between the characteristic time of collapse and the characteristic time of heat transfer.
Adiabatic collapse. The average pressure and the temperature in the vapour cavity increase during the collapse (see Fig 11), locally reaching a max P � 40 MPa and T � 10000 K. Those values, despite some difference in the operating conditions, are comparable to those measured by Obreschkow et al. [72].
The pressure and temperature field distribute differently in the cavity: 1. In the first stage of the collapse, the interaction between gas and fluid results in a rapid increment of pressure at the cavity interface (Fig 12a) generating a shock wave inside the cavity. Later, because of the combined effect of cavity compression and shock wave propagation, the pressure increases in the centre of bubble (Fig 12b) becoming almost uniform at the collapse.
2. Similarly to the pressure, the temperature increases at the cavity interface during the first stages of the collapse (see Fig 12c). The absence of heat diffusion mostly affects the final stage of the collapse: the internal energy does not diffuse and heat is confined and accumulated at the cavity interface (Fig 12d).
(Heat) Diffusive collapse. Our results show that time of the collapse is not significantly affected by the presence of heat transfer in the model. However, the pressure and temperature peak inside the bubble decreases with respect to the adiabatic case, see Fig 13. Also the pressure and temperature fields inside the cavity change: 1. In the first phase of the collapse, the pressure field in the cavity remains uniform (Fig 14a).
Approaching the final stage of the collapse, the pressure is slightly lower at the cavity interface, Fig 14b, because of the presence of the diffusive heat transfer mechanism. In fact, the pressure of an ideal gas is function of both the density and the internal energy, as shown by Eq 17. The presence of a high temperature gradient at the cavity interface greatly reduces the internal energy of the particles in that area.
2. Similarly to the adiabatic case, in the first stage the interaction between gas and the fluid tend to increment the temperature near the cavity interface, see Fig 14c. However, in this case the heat generated at the interface diffuses both towards the centre of the cavity and into the liquid. This leads to higher temperatures at the centre than at the interface (see Fig 14d).
By comparing Fig 14 with Fig 12, it is clear that, despite the total collapsing time is almost the same, the heat exchange mechanism has important consequences on both temperature and pressure in the cavity. Our results, therefore, show that adiabatic conditions, despite being often used in the literature [30,31,[73][74][75], are not always realistic. The collapse generates great amount of heat and high temperatures are reached in the cavity. Despite the small timescale of the process, temperature in the cavity rapidly grows and the heat transfer from the cavity interface to the surroundings cannot be neglected.
Effect of thermal diffusivity on the temperature peak. In the previous section, we calculated a specific case, where the liquid is water, the gas is water vapour, ΔP = P L − p G = 5MPa and R 0 = 100μm. In this section, we study how different parameters and initial conditions would affect the temperature rise T/T 0 in the cavity. In order to simplify the study, we perform

PLOS ONE
a dimensional analysis of our system to reduce the number of significant parameters. Assuming that the collapse depends on ΔP, ρ L , α L,G = (α L + α G )/2 and R 0 , and using the Buckingham π theorem, it is possible to determine that the system depends on two fundamental dimensionless groups

PLOS ONE
The dimensionless group P 1 can also been seen as the ration between the characteristic heat diffusion time scale of the process, R 2 /α L,G , and the Rayleigh collapsing time, R ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi r L =DP p : When t c � τ d , the collapse is faster than the characteristic time of heat transfer, the heat generated is trapped in the cavity, and the process can be considered adiabatic. When t c � τ d , the characteristic time of heat transfer is smaller that the collapsing time and the process can be considered isotherm.
In Fig 15 T

Conclusion
This work proposes the first SPH model of a collapsing cavity filled with non condensable gas coupled with the heat transfer mechanism.
The aim of the work is to understand the role of diffusive heat transfer during the Rayleigh collapse. This was achieved by introducing the dimensionless group, P 1 . P 1 defined as ratio between the characteristic time of collapse and the characteristic time of thermal diffusion.
In Fig 15 five regions are identified. For each of these regions the temperature field of the gas-liquid system distributes differently:

PLOS ONE
• I region (0 < P 1 < 6.5 � 10 −3 ): in this region, both the gas and liquid behave isothermally. As a result of this, the liquid drains all the energy developed by the gas.
• II region (6.5 � 10 −3 < P 1 < 10): in this region, the gas behaves isothermally while the liquid shows a temperature profile. The energy rapidly diffuses inside the cavity, flattening the temperature profile in the cavity.
• III region (10 < P 1 < 1 � 10 5 ): in this region, neither the gas nor the cavity are adiabatic. This scenario is described in (heat) diffusive collapse section and by Fig 14. • IV region (1 � 10 5 < P 1 < 5 � 10 7 ): in this region, the gas in the cavity behaves adiabatically, but the liquid does not. Therefore, part of the energy generated at the interface is transferred to the liquid phase.
• V region (P 1 > 5 � 10 7 ): in this region both the gas in the cavity and the liquid behave adiabatically. All the energy generated by the collapse is trapped in the cavity and the temperature increment is concentrated in the cavity interface (see Fig 12).
In brief, this analysis shows that for P 1 > 5 � 10 7 the collapse can be considered adiabatic. At smaller P 1 , the heat generated at the cavity interface is taken away by the liquid phase, or diffuses in the cavity, homogenising the temperature field. When P 1 < 6.5 � 10 −3 , all the heat generated in the collapse is drained by the liquid and the collapse can be considered isotherm.
This shows that, despite the short timescale, the presence of the heat transfers mechanism leads to a temperature peak drop of around 50% compared to the adiabatic case. This suggest that the adiabatic assumption for the Rayleigh collapse could leads to a not reliable pressure and temperature peak estimation and its use should be properly justified.