Numerical simulation of cavitation-bubble expansion and collapse inside a bottle subjected to impact on its topside

Abstract

In this study, we simulate the acceleration-induced cavity bubble expansion and its non-spherical collapse in a drastically varying pressure field caused by mechanical impact, and we qualitatively compare the results with those for similar ζ number (Obreschkow et al., 2011) as recorded by a high-speed camera (Taylor, 2012). The gas–liquid interface is tracked by using the volume of fluid method. The nonlinear compressibility effect of the liquid phase is considered by applying the Tait equation. The simulation of the impact load is accomplished by converting the acceleration curve obtained from the experiment (Kang & Raphael, 2018) into the velocity curve of the bottom boundary movement through integration and subsequent using of dynamic mesh model. Numerical results show that cavitation bubble nucleuse expand rapidly in the drastically decreasing pressure field caused by the impact and a pressure wave is emitted and propagates in the liquid during the subsequent collapse of the bubble. When the pressure wave is reflected on the wall, it causes a high peak pressure at the wall surface, which may underlie the rupture of the bottle wall in the experiment (Daily & Pendlebury, 2014).

KEYWORDS:

• High acceleration
• cavitation
• impact load
• dynamic mesh
• bubble collapse
• pressure wave

1. Introduction

Striking the top of a bottle filled with water can cause the water at the bottom of the bottle to experience cavitation, which can subsequently cause the bottle wall to crack (Daily & Pendlebury, 2014). This cavitation induced by acceleration is different from that induced by the fluid momentum, which is governed by the conventional cavitation number. Therefore, Pan et al. proposed a new cavitation number, $Ca$, and experimentally verified that it can be applied as a general criterion for the onset of acceleration-induced cavitation (Fatjo, 2016; Pan et al., 2017). As well, many scholars have studied the dynamics of bubble impacted by a shock wave, and these works have demonstrated the jetting mechanism of the bubble if exposed to a shock wave (Hawker & Ventikos, 2012; Ohl and Ikink, 2003; Ohl & Ohl, 2013; Seung & Kwak, 2017). Although the criteria for the onset of cavitation under impact have been clearly established and the dynamics of bubble impacted by a shock wave have been described, the dynamic of the gas–liquid interaction in a bottle under mechanical impact is still largely unknown. Accurately describing and predicting the spatio-temporal dynamics of bubbles under shock-loading conditions are essential for many application scenarios such as pipeline flow and biological tissue damage from blast or impact exposures (Hong et al., 2016; Kang & Raphael, 2018; Rosenfeld et al., 2013). Therefore, it is particularly important to understand the dynamic mechanism underlying gas–liquid interactions in a bottle under an impact load.

In numerical research, many methods such as the boundary element method, level-set method, Eulerian finite element method and volume of fluid (VOF) method have been used to describe the dynamic behavior of the bubble (Akhatov et al., 2001; Lauer et al., 2012; Liu et al., 2018; Liu et al., 2020; Luo, 2008; Méndez & González-Cinca, 2011; Peng et al., 2019; Tang et al., 2020). Among these, the VOF method was used by Koukouvinis to study the aspheric collapse of bubble subject to gravity (Koukouvinis et al., 2016a, 2016b). Another notable example is the work of Osterman et al., which used the VOF method to simulate the near-wall bubble collapse in an ultrasonic field(Osterman et al., 2009). In this study, taking into account the characteristics of acceleration-induced cavitation, we used the VOF and dynamic mesh method to study the expansion and collapse behavior of acceleration-induced cavitation under mechanical impact. We also examine the subsequent formation and propagation of pressure waves.

In a previous experimental study, a non-dimensional scaling law was proposed by a team from École polytechnique fédérale de Lausanne (EPFL) for characterizing the formation of jets from bubble collapse, in which a dimensionless number $\zeta$ was used to characterize the collapse asymmetry (Obreschkow et al., 2011): (1) $\zeta =\frac{|\mathrm{\nabla }p|R_{max}}{p_{\mathrm{\infty }}-p_v}$(1) Here, $\mathrm{\nabla }p$ denotes the pressure gradient due to acceleration; $R_{max}$ is the maximum bubble radius; $p_{\mathrm{\infty }}$ is the pressure at the same horizontal position within the bubble; and $p_v$ is the vapor pressure. When $\zeta \ge 0.0004$, the jet pierces the bubble. In this case, the maximum instantaneous acceleration is 500 g ($g=9.81m/s^2$), corresponding to a instantaneous ζ value of about 0.05. According to this rule, we can predict that significant non-spherical collapse and jetting will occur in this example (Obreschkow et al., 2013).

2. Numerical model

The commercial flow solver Fluent17.1 is used for the numerical simulation, wherein the VOF method is applied to accurately track the gas–liquid interface. Surface tension effects are considered and the surface tension coefficient is set to $0.072N/m$, although its influence during jet formation is negligible as the Weber number in this process is ∼ 3000. The mass transfer between the two phases is ignored. The equation solved, based on the Navier-Stokes equation in compressible viscous form, are as follows:

Continuity equation: (2) $\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (\rho \overrightarrow{v})=0$(2) where $\rho$ is the density of the fluid and $\overrightarrow{v}$ is the velocity vector of the flow field.

Momentum equation: (3) $\frac{\mathrm{\partial }}{\mathrm{\partial }t}(\rho \overrightarrow{v})+\mathrm{\nabla }\cdot (\rho \overrightarrow{v}\overrightarrow{v})=-\mathrm{\nabla }p+\mathrm{\nabla }\cdot (\mathit{\tau })+\rho \overrightarrow{g}+\overrightarrow{f}$(3) where $p$ is the pressure, $\overrightarrow{g}$ is the gravity vector, $\overrightarrow{f}$ are body forces, and $\mathit{\tau }$ is the viscous stress tensor, defined as follows: (4) $\mathit{\tau }=\mu [(\mathrm{\nabla }\overrightarrow{v}+\mathrm{\nabla }{\overrightarrow{v}}^T)]-\lambda \mathrm{\nabla }\cdot \overrightarrow{v}I$(4) where $\text{μ}$ is the dynamic viscosity of the fluid, $\text{I}$ is the identity matrix affected by bulk expansion, and $\lambda$ is the bulk viscosity of the fluid, usually set to $-2/3\text{μ}$. The Reynolds number of the flow varies from less than 1000 during the bubble expansion to nearly 40,000 in a short time period of bubble collapse. The flow regime was in a transitional state of laminar to mildly turbulent flow, with laminar flow for most of the time. And the mesh used in the simulation is dense enough to be considered as a direct numerical simulation (DNS) of the flow process, so no turbulence model is used.

Energy conservation equation: (5) $\begin{array}{rl}& \frac{\mathrm{\partial }}{\mathrm{\partial }t}(\rho E)+\mathrm{\nabla }\cdot (\overrightarrow{v}(\rho E+p))\\ & =\mathrm{\nabla }\cdot \left(k\mathrm{\nabla }T-\sum _j{h}_j{J}_j+(\mathit{\tau }\cdot \overrightarrow{v})\right)+S_h\end{array}$(5) where $E=h-\frac{p}{\rho }+\frac{v^2}{2}$ denotes the total energy per unit mass and the first three terms to the right of the equation represent the energy transfer due to conduction, species diffusion and viscous dissipation, respectively. The last term represents the volumetric heat source and according to the two-phase model chosen for this study, it's zero.

Volume fraction equation: (6) $\frac{\mathrm{\partial }{\alpha }_i{\rho }_i}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot ({\alpha }_i{\rho }_i\overrightarrow{v})=0\text{, }i=1,2$(6) where ${\alpha }_i$ is the volume fraction of the fluid and the subscript $i$ indicates the primary or secondary phase. The interface between the two phases can be determined by solving the equation for the second phase.

Both phases are considered compressible, where the ideal gas model is used for the gas and the following Tait equation is used for the liquid: (7) $p=\frac{{\rho }_0{c_0}^2}{n_l}({\left(\frac{\rho }{{\rho }_0}\right)}^{n_l}-1)+p_0$(7) where ${\rho }_0$ is the density of water, equal to $998.2kg/m^3$, $c_0$ is the speed of sound, equal to $1450m/s$, at the reference state $p_0=3540Pa$ and the exponent $n_l$ is set to 7.15. The Tait equation accounts for a nonlinear compressible effect, which is a prerequisite for tracking pressure waves.

The finite volume method was used to solve the Naver-Stokes equations for the transient state on a 2D grid, with the SIMPLE algorithm for pressure-velocity coupling and the implicit scheme for time discretization. A second-order upwind scheme was used for the discretization of density and momentum. The convergence criteria were ${10}^{-5}$ for the continuity and momentum equations, and ${10}^{-9}$ for the energy equation. Time stepping was done with an adaptive method, controlling the Courant number to be less than 0.2, with a variable time step of ${10}^{-9}\mathrm{s}$ to ${10}^{-7}\mathrm{s}$. Hawker and Ventikos examined the three-dimensional effect in detail in their study of the interaction of a strong shock wave with a gas bubble in a liquid medium (Hawker & Ventikos, 2012). The results show that, apart from small-scale details in the collapse phase, there is a lack of significant three-dimensional structure during the bubble evolution, and in addition the collapse period, maximum jet velocity and maximum pressure differ by about 15-25% between the two computation. Since 2D simulation significantly reduces the computational cost as compared to the full 3D simulation, the 2D computational model was chosen here.

In addition to the solution of the Navier-Stokes equations, a widely used tool for predicting bubble expansion and collapse in a varied pressure field is the Keller-Miksis equation, which takes the following form: (8) $\begin{array}{rl}& (1-\frac{U_b}{C})R_b{\dot{U}}_b+\frac{3}{2}{U}_b^2(1-\frac{U_b}{3C})\\ & =(1+\frac{U_b}{C})\frac{P_l-p_0}{{\rho }_l}+\frac{R_b}{{\rho }_{l}C}\frac{d(P_l-p_0)}{dt}\end{array}$(8) Where $R_b$ is the bubble radius, $U_b=d{R}_b/dt$ means the velocity of the bubble wall, $C$ is the speed of sound in the liquid, $p_0$ and ${\rho }_l$ are the pressure and density of the ambient liquid respectively, $P_l$ represents the liquid pressure at the bubble wall. The gas inside the bubble obeys the van der Waals law, then $P_l$ is given as follows: (9) $\begin{array}{rl}{P}_l& =(p_0-p_v+\frac{2\sigma }{R_0}){\left(\frac{R_0^3-b{R}_0^3}{R_b^3-b{R}_0^3}\right)}^n\\ & -\frac{2\sigma }{R_b}-\frac{4\mu U_b}{R_b}+p_v\end{array}$(9) where $p_v$ is the vapor pressure and ignored in this study, $R_0$ is the radius of the bubble at equilibrium state, $\mu$ and $\sigma$ are the viscosity and surface tension of the liquid. $n$ is the polytropic exponent and $b$ is the van der Waals parameter.

3. Validation of the methodology

A canonical case of an isolated spherical bubble collapsing in the free-field has been used for the validation of the methodology, as shown in Figure 1. The initial radius of the bubble is $10\mu m$, the internal pressure is $6900Pa$ and the ambient pressure is $1Bar$. The radius of the spherical calculation domain is 1000 $\mu m$, 100 times the initial radius of the bubble, to minimize the influence of the boundary on the evolution of the bubble. Comparison of the Keller-Miksis equation with the numerical results is shown in the Figure 2. The results of fine mesh were in good agreement with the analytical solution, indicating that the adopted numerical model could be used to simulate the dynamic process of a single bubble. In addition, the mesh adopted in the following case study kept the same level as the fine mesh used in the validation to ensure the accuracy of the result.

Figure 1. Computational domain for the validation.

Figure 2. Comparison of the results of Keller-Miksis equation with numerical calculations.

4. Computational mesh and boundary condition

As shown in Figure 3, the computational region is 80 mm long and 40 mm wide. The mapped-type structured mesh and the regional mesh adaption technique are used to track the bubble interface, particularly at the moment of collapse. The total number of grids is about 680,000. Particularly in the bubble-collapse region, the minimum grid size is reduced from $10\mu m$ to $0.625\mu m$ by means of the regional mesh adaption technique.

Figure 3. The 2d computational domain and the mapped-type structured mesh with refinement in the bubble-collapse region.

The accurate tracking and 'capture' of impact information is the key to the successful simulation of impact loads. In this regard, Kang et al. designed and used a drop-tower-based integrated system, which can suitably capture acceleration with time under mechanical impact (Kang et al., 2017b). In this study, we apply and integrate the curve fitting of the acceleration data obtained from the experiment to obtain the velocity change curve of the bottle under impact. This procedure makes it possible to use the dynamic mesh model to simulate impact loads in numerical simulations. This use of moving boundaries to obtain changing pressure fields is also common in bubble simulations under the influence of acoustic fields (Osterman et al., 2009; Boyd & Becker, 2019).

To facilitate the curve integral operation, we selected a set of sinusoidal functions to fit the curve of acceleration, which can be expressed as follows: (10) $\mathrm{A}(t)=\sum _{n=1}^6{a}_n\ast \sin (b_n\ast t+c_n),(0<t<0.0012)$(10) In the dynamic mesh model set-up, the numerical simulation of the dramatic pressure field caused by the impact was achieved by controlling the velocity of motion over the entire area, and the velocity function can be obtained by integrating the above formula. It can be expressed as follows: (11) $\mathrm{V}(t)=-\sum _{n=1}^6\frac{a_n}{b_n}\ast \cos (b_n\ast t+c_n)+d,(0<t<0.0012)$(11)

Here, $a_n$, $b_n$, $c_n(1\le n\le 6)$, $d$ are listed in Table 1. And the curves of acceleration and velocity are shown in Figure 4.

Figure 4. Left: Acceleration fitting curve after the bottle is impacted. Right: Velocity curve obtained by integrating acceleration.

Table 1. Coefficients of the fitting curve function.

i	ai	bi	ci	d
1	2394	5137	2.367	−0.5742
2	2066	1575	−1.785
3	777.4	9027	1.046
4	276.4	1.934e+04	−2.241
5	518.4	2.322e+04	−1.926
6	455	2.364e+04	0.9507

In one calculation example, the computational region corresponds to only water at the beginning, and the dynamic boundary conditions are used to obtain the pressure field, which is used as a comparative calculation example when no cavitation occurs. In another example, with the application of the moving boundary condition, when the pressure at the pre-set cavitation point $P$ ($0.02m$, $0.003m$) is less than the vapor pressure ($P_v=3540Pa$ at temperature $300K$), a vapor bubble nucleus with a radius of $R_0=100\mu m$ is introduced. The case of a bubble nucleus height of $0.003m$ was chosen as the typical case to be presented with details in this paper, since the maximum pressure drop occurs at the bottom of the bottle when the impact load is applied and the scenario when the bubble expands to its maximum was evaluated to make sure that it does not come into contact with the bottom.

5. Result analysis

When compared with the pressure field without cavitation, the presence of a cavitation bubble significantly reduces the amplitude of the pressure drop and delays the time for the pressure to return to the initial state, as can be seen from Figure 5. In the case of impact loading, the presence of a small cavitation nucleus ($100\mu m$) will significantly alter the pressure changes across the entire flow field and play an important role in buffering the pressure drop.

Figure 5. Comparison of the pressure at position (0.01, 0.003 m) for different cases with and without cavitation onset.

The delay effect of the pressure recovery increases the time required for bubble collapse, which is consistent with the observation that a single bubble is affected by the surrounding bubble cloud and the collapse time is greatly delayed (Taylor, 2012; Daily & Pendlebury, 2014). This situation also affects the comparison and verification of simulation results with the experiments. Because the ambient pressure field of the bubble observed in the experiment is affected by the bubble clouds, resulting in a deviation in the evolution time between the simulation and the experiment. Therefore, here, we only qualitatively compare the evolution mode of the bubble. That is, the comparison of bubble shapes during the bubble collapse phase with a dimentionless time, which is obtained by dividing the evolution time by the total collapse time, as shown in Figure 6. In a high-speed video recording of the bubble collapse (Taylor, 2012), we selected a relatively independent bubble at a certain distance from the bubble-cloud area at the bottom of the bottle, and the bubble contour obtained by computational fluid dynamics (CFD) simulation was compared with the image captured in the experiment. The results showed that the simulation and experiment exhibit a high degree of consistency as regards the topological changes of the shape, which proves the reliability of the numerical calculation.

Figure 6. Shape comparison of simulation results with experimental images during the bubble collapse phase. The total bubble collapse times in the experiment and simulation are 0.5, 0.1 ms respectively.

When the bottle is impacted and suddenly accelerates, the water in the bottle is 'stretched' due to inertia, resulting in a drastic pressure drop. When the bottom pressure of the bottle is less than the vapor pressure, vapor-bubble nucleuses appear (Schnerr & Sauer, 2001; Fatjo, 2016; Kang et al., 2017a). As the acceleration of the liquid further increases, the pressure continues to decrease, and the pressure difference between the bubble and the surrounding liquid becomes the main driving force underlying bubble growth and collapse. In the initial stage of bubble expansion, the bubble remains spherical. In the late expansion stage, the shape of the bubble deviates from the spherical shape, and the bubble is slightly flattened in the lateral direction owing to the pressure difference as shown in Figure 7. Finally, the expansion of the bubble causes a significant drop in pressure within the bubble. Simultaneously, the ambient pressure continues to rise, which eventually slows down the bubble expansion process and causes its subsequent collapse.

Figure 7. Contour evolution diagram of the bubble driven by the changing ambient pressure at position (0.01 m,0.003 m).

When the bubble radius reaches the maximum value, it begins to collapse. In the initial stage of collapse, along the direction of acceleration, there are minor differences in the rate of collapse between the upper and lower bubble walls. In addition to the non-uniform distribution of pressure around the bubbles caused by the rigid wall, the ultra-high acceleration caused by the impact also results in a gradient distribution of the pressure field. And due to the higher pressure acting on the top of the bubble, the collapse velocity of the upper wall is greater than that of the lower wall, and thus, a recessed area is formed here. As the bubble collapses further, this deviation will cause a phenomenon of momentum focusing and form a positive feedback mechanism, resulting in even greater velocity at the top of the bubble (Koukouvinis et al., 2016a). Eventually, the bubble approaches its minimum radius, and a jet in the opposite direction to the pressure gradient forms at the top of the bubble.

In Figure 8, a jet flow develops gradually from a small concave at the top of the bubble and penetrates the bubble from the inside. As the jet impacts the opposite bubble wall, a high-pressure zone is generated, and there is a ring-shaped high-pressure expansion wave above and below the zone, where the first pressure wave, i.e. a water hammer shock, is formed.

Figure 8. Scenario of water hammer shock formed by jet impinging the wall. The black line indicates the interface between gas and liquid (a) Jet flow inside of the bubble, (b) The high-pressure zone generated by the jet piercing the bubble wall, (c) Pressure and velocity field after the jet penetration of the bubble wall.

During impact, the blunt rounded front of the jet contacts the lower bubble wall creating a very sharp lower rim of the bubble. As the jet continues to impact the lower wall, the gas trapped at the lower rim of the bubble is further squeezed to the sides of the remaining bubble. The 'expelled-gas shock' appears inside the lower side of the bubble (the red section within the contour line in Figure 8(c)). A similar phenomenon also occurs in the laser-induced bubble collapse stage near a rigid wall (Lechner et al., 2017). At this time, the jet radius is $230\mu m$, larger than the initial bubble radius, and the velocity is $23m/s$. With the widening of the liquid jet and the propagating of the high-pressure expansion wave, the impact process of the jet comes to an end and transforms into the collapse phase of the torus bubble (3D view), which is accompanied by the formation of secondary jets and the emission of pressure waves.

In Figure 9, at the lower rim of the torus bubble, a small, point-like high pressure zone drives a tiny annular jet from its lower edge into the bubble. Subsequently, this nanojet enters the bubble at an oblique upward angle from the central axis. When compared with the main jet, in addition to being significantly smaller in size, the nanojet has another different feature, that is, the main jet is an axial jet, while it is an annular jet. This annular jet soon strikes the outer wall of the torus bubble, and another pressure wave is emitted. The laser-induced bubble collapse near a rigid wall also exhibits similar features (Lechner et al., 2017). The collapse of the torus bubble progresses from the lower rim upwards and is accompanied by the overall contraction of the bubble.

Figure 9. Collapse stages of the torus bubbles at different times (a) A nanojet is driven by a small, point-like high pressure zone, (b) A nanojet enters into the bubble, (c) Emission of the secondary pressure wave.

When the torus bubble finally collapses, a ring of high pressure is accompanied. When the pressure wave propagates to the axis of symmetry, the local pressure reaches a peak due to the superposition of the oppositely propagating impact wave. Subsequently, the pressure along the axial direction remains high value, propagating upward to the free surface and downward to the bottom wall, as shown in Figure 10.

Figure 10. Superposition of the pressure wave.

In Figure 11, the downward propagating pressure wave is reflected on the flat surface of the lower wall and propagates in the opposite direction while geometrically maintains the same shape as before, that is, a pressure wave that spreads around in a folded-spherical way. Subsequently, it strikes the side wall surface and is reflected again.

Figure 11. Propagation and reflection of the pressure wave inside the bottle (a) Reflection of the pressure wave at the lower wall, (b) Propagation of pressure wave in the liquid, (c) Reflection of the pressure wave at the side wall.

We focus on the pressure at the bottom wall, especially the pressure change at the point of crossover with the axis of symmetry, where the pressure is the highest because it is located at the 'apex' of the torus bubble. In Figure 12, the first peak of the blue line is the result of the collapse of the torus bubble and the pressure wave propagating to the lower wall. At this time, the pressure peak exceeds 100 bars. The second peak of the blue line is the result of the pressure wave passing the midpoint of the bottom surface after being reflected off the side-wall surface, and it is lower in magnitude than the first peak. The first peak of the green line is the result of the first pressure wave propagating to the side wall, with a peak pressure of 46 bars, which is less than the value of the peak at the midpoint of the bottom wall. Subsequently, the pressure wave reflects back and forth off the left and right walls, and the peak value is gradually attenuated.

Figure 12. Pressure curve at the lower-wall and side-wall surfaces.

6. Discussion and conclusion

In the case considered here, when the bottle is subjected to a mechanical impact, the maximum acceleration reaches $500g$, corresponding to an instantaneous $\zeta$ value of about 0.05, which provides a rare study case of bubble collapse in an extremely highly accelerated field. Based on the numerical results, we quantitatively describe the dynamic response of a cavitation bubble in a drastically changing pressure field when subjected to an impact load, and we analyze the formation and propagation of the two pressure waves resulting from the bubble collapse.

We compared the different cases with and without cavitation, and our results show that the presence of cavitation significantly alters the pressure changes across the entire flow field, reduces the amplitude of the pressure drop, delays the time for the pressure to return to the initial state, and overall plays a significant role in buffering the pressure change. This is because the impact load transforms the liquid in the bottle from a pure liquid state to a gas–liquid mixed state. When compared with pure liquid, this mixed gas–liquid state plays a role in buffering the rigidity of the liquid and increasing the overall elasticity (stretchability) of the fluid.

Simultaneously, the delay effect of pressure recovery also significantly prolongs the time of the bubble collapse stage, which leads to a deviation between simulation and experiment. Because the observed bubble lies near a bubble cloud, its ambient pressure is affected by the bubble cloud, and thus, there are deviations in the evolution time between simulation and experiment, and only a qualitative comparison can be performed in terms of the evolution mode of the bubble.

From the perspective of energy, along with the expansion stage of the bubble, the mechanical energy of the impact is transformed into the potential energy of the bubble. When the bubble collapses, most of the bubble energy radiates outwards in the form of pressure waves and is finally transferred to the nearby wall surface. This mechanism may underlie the rupture of the bottle wall in the experiment.

Acknowledgements

The authors would like to acknowledge the financial grants from the National Project of China (No. 6140206040301).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Project of China [Grant Number 6140206040301].