Modeling Acoustic Cavitation Using a Pressure-Based Algorithm for Polytropic Fluids

Abstract

A fully coupled pressure-based algorithm and finite-volume framework for the simulation of the acoustic cavitation of bubbles in polytropic gas-liquid systems is proposed. The algorithm is based on a conservative finite-volume discretization with collocated variable arrangement, in which the discretized governing equations are solved in a single linear system of equations for pressure and velocity. Density is described by the polytropic Noble-Abel stiffened-gas model and the interface between the interacting bulk phases is captured by a state-of-the-art algebraic Volume-of-Fluid (VOF) method. The new numerical algorithm is validated using representative test-cases of the interaction of acoustic waves with the gas-liquid interface as well as pressure-driven bubble dynamics in infinite and confined domains, showing excellent agreement of the results obtained with the proposed algorithm compared to linear acoustic theory, the Gilmore model and high-fidelity experiments.

1. Introduction

Acoustic cavitation describes the pressure-driven behavior of bubbles in a liquid environment, where an externally introduced change in pressure, e.g., by an ultrasound transducer or through laser-induced optical breakdown, causes a radial bubble motion that can range from moderate bubble oscillations to a strong inertial bubble collapse that emits shock waves in the liquid, concentrating large amounts of acoustic energy [1]. In many engineering applications, ranging from fast-running ship propellers to artificial heart valves, hydrodynamic cavitation, where the pressure difference is induced by an acceleration of the flow, is associated with considerable negative effects on its environment and is, therefore, undesirable. Acoustic cavitation, however, was found to be useful for a large number of engineering applications, including ultrasonic cleaning [2], for ultrasonic contrast agents in medical imaging [3], to actuated microbubbles that act as micropumps [4], and to facilitate self-organization of microbubbles into crystals and arrays [5].

Since Lord Rayleigh first derived a model for the pressure-driven collapse of an empty cavity in an incompressible liquid in 1917 [6], substantial research efforts have been dedicated to understanding and mathematically describing the pressure-driven behavior of bubbles [1,7]. Plesset [8] extended the model of Rayleigh to account for the gas phase inside the bubble, which led to the famous Rayleigh-Plesset (RP) equation, a second-order differential equation, which has since been extended to include liquid compressibility [9,10], among other mechanisms. The RP equation and its subsequent extensions have been the basis for most subsequent research on cavitation. However, RP equations are based on strongly simplifying assumptions: the bubble is spherical and the gas density is assumed to be much smaller than the liquid density. Alleviating these limitations requires computation of the fully resolved flow field and interface dynamics of the gas-liquid system, by solving the full set of governing conservation laws describing the behavior of two-phase flows using appropriate numerical methods, such as finite-volume or finite-element methods. Such numerical methods have been recently applied to study a broad range of cavitation phenomena, including laser-induced cavitation bubbles [11,12], the bubble-collapse-driven penetration of a biomaterial surrogate [13], shock-driven bubble collapse [14,15] or the dynamics of cavitation bubbles in viscoelastic media [16].

Describing the interacting fluids by a polytropic equation of state is a popular choice when modeling the pressure-driven behavior of bubbles [7,12]. In a polytropic fluid, the density $\rho$ is a function of only the pressure p, as defined by the polytropic relationship, given for an ideal gas as $\rho =K{p}^{1/\kappa }$, where $\kappa$ is the polytropic exponent and $K>0$ is a fluid-dependent constant. From a mathematical and numerical viewpoint, the most appealing feature of flows of polytropic fluids is that pressure does not depend on temperature and the energy equation becomes redundant, simplifying the mathematical analysis of these flows as well as the design and implementation of the associated numerical methods. Describing the gas inside bubbles as a polytropic fluid is one of the founding assumptions of the RP equation [17] and using polytropic fluid models has been demonstrated to be suitable even for the prediction of complex bubble behavior, such as wall-bounded cavitation [12].

Pressure-based algorithms, in which the density is described by a suitable equation of state and the continuity equation serves as an equation for pressure, are particularly suited to model acoustic cavitation, due to the broad range of Mach numbers occurring in cavitating gas-liquid systems. To this end, an important aspect of pressure-based algorithms is that changes in pressure are finite in all Mach number regimes and that the vanishing density differences at low Mach numbers do not present a problem [18]. Pressure-based algorithms for compressible and incompressible flows are traditionally founded on a predictor-corrector procedure, e.g., projection methods [19,20] or the SIMPLE method [21]. Previously proposed pressure-based finite-volume algorithms using polytropic fluid models dedicated to the prediction of bubble dynamics and cavitation are also founded on such segregated algorithms [12,22,23]. Over the past decade, research efforts dedicated to pressure-based algorithms have increasingly focused on fully coupled methods [18,24,25,26,27,28,29], whereby the discretized governing equations are solved simultaneously in a single system of equations, which have been shown to yield performance benefits for incompressible flows [24] and an improved stability for two-phase flows [26]. Furthermore, fully coupled algorithms do not require any form of underrelaxation to reach a converged solution, even for compressible flows with strong shocks [29], contrary to segregated pressure-based algorithms. Recent work by Denner and co-workers [15,30] has demonstrated the utility of fully coupled pressure-based algorithms for interfacial flows at all Mach numbers, including bubble dynamics and collapse. Despite the demonstrated versatility, robustness and increasing maturity of fully coupled pressure-based algorithms, such an algorithm has not yet been proposed for two-phase flows of polytropic fluids.

In this article, a novel fully coupled pressure-based algorithm for the prediction of acoustic cavitation is proposed, based on a conservative finite-volume framework and a Volume-of-Fluid (VOF) method to represent the interacting bulk phases. The discretized continuity and momentum equations are solved simultaneously in a single linear system of equations for pressure and velocity, with density defined based on pressure by a polytropic fluid model. Each term of the linearized and discretized governing equations makes an implicit contribution to pressure and/or velocity, leading to a robust pressure-velocity coupling. The proposed numerical framework is validated by representative test-cases, including the propagation of acoustic waves, as well as the Rayleigh collapse and the wall-bounded cavitation of a bubble.

The governing equations are presented in Section 2, the polytropic closure model is presented in Section 3, the numerical framework is described in Section 4 and the interface treatment is presented in Section 5. The results of the various test-cases considered to validate the numerical framework are presented in Section 6. The article is concluded and summarized in Section 7.

2. Governing Equations

The conservation laws governing the flow of a polytropic fluid are the conservation of mass

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho u_i}{\partial x_i}=0$$

(1)

and the conservation of momentum

$$\frac{\partial \rho u_j}{\partial t}+\frac{\partial \rho u_i{u}_j}{\partial x_i}=-\frac{\partial p}{\partial x_j}+\frac{\partial {\tau }_{ji}}{\partial x_i},$$

(2)

where t is the time, $\mathit{u}$ is the velocity vector, p is the pressure and $\rho$ is the density. The stress tensor $\tau$ for the considered Newtonian fluids is given as

$${\tau }_{ji}=\mu \left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right)-\frac{2}{3}\mu \frac{\partial u_k}{\partial x_k}{\delta }_{ij},$$

(3)

where $\mu$ is the dynamic viscosity.

The Volume-of-Fluid (VOF) method [31] is applied to distinguish the interacting bulk phases and model the motion of the fluid interface. Both bulk phases are represented by an indicator function $\zeta \left(x\right)$, with

$$\zeta \left(x\right)=\left\{\begin{array}{cc}0\hfill & \mathrm{if}x\in {\mathsf{\Omega }}_{\mathrm{a}}\hfill \\ 1\hfill & \mathrm{if}x\in {\mathsf{\Omega }}_{\mathrm{b}}\hfill \end{array}\right.$$

(4)

where $\mathsf{\Omega }={\mathsf{\Omega }}_{\mathrm{a}}\cup {\mathsf{\Omega }}_{\mathrm{b}}$ is the computational domain, where ${\mathsf{\Omega }}_{\mathrm{a}}$ and ${\mathsf{\Omega }}_{\mathrm{b}}$ are the subdomains occupied by fluids a and b, respectively. Neglecting mass transfer between the bulk phases, the fluid interface propagates with the flow and, therefore, the material derivative of $\zeta$ is

$$\frac{D\zeta }{Dt}=\frac{\partial \zeta }{\partial t}+u_i\frac{\partial \zeta }{\partial x_i}=0.$$

(5)

In the interest of clarity but without loss of generality, surface tension and gravity are neglected in the following, to focus on the novel aspects of the proposed algorithm. However, the extension of the proposed algorithm to include surface tension and gravity is straightforward, following the work of Denner and van Wachem [26].

3. Polytropic Closure

The governing Equations (1) and (2) are closed by defining the density based on pressure, using the Noble–Abel stiffened-gas (NASG) model [32]. The NASG model combines the Noble–Abel gas model (also often called co-volume gas model) [33] and the stiffened-gas model [34] to obtain a simple gas model that accounts for molecular repulsion and attraction. The thermal equation of state of the NASG model is defined as [32]

$$p(v,T)=(\gamma -1)\frac{c_{v}T}{v-b}-\mathsf{\Pi },$$

(6)

where $\gamma =c_p/c_v$ is the heat capacity ratio, $c_v$ is the specific isochoric heat capacity, $c_p$ is the specific isobaric heat capacity, $v=1/\rho$ is the specific volume, b is the co-volume and $\mathsf{\Pi }$ is a pressure constant that accounts for the attraction between molecules. Based on the isentropic relationship for such an NASG fluid [32], a polytropic fluid model can be derived as

$$\frac{\rho }{1-\rho b}=K{(p+\mathsf{\Pi })}^{\mathsf{\Gamma }},$$

(7)

where $\mathsf{\Gamma }=1/\kappa$ and $\kappa$ is the polytropic exponent, from which the density follows as

$$\rho =\frac{K{(p+\mathsf{\Pi })}^{\mathsf{\Gamma }}}{1+bK{(p+\mathsf{\Pi })}^{\mathsf{\Gamma }}}.$$

(8)

The polytropic constant K is defined based on the predefined reference density ${\rho }_0$ and the predefined reference pressure $p_0$ as

$$K=\frac{{\rho }_0}{{(p_0+\mathsf{\Pi })}^{\mathsf{\Gamma }}(1-b{\rho }_0)}.$$

(9)

The polytropic NASG model reduces to the polytropic description of an ideal-gas model for $\mathsf{\Pi }=0$, $b=0$, of a Noble-Abel gas model for $\mathsf{\Pi }=0$, $b>0$, and of a stiffened-gas model for $\mathsf{\Pi }>0$ and $b=0$.

4. Numerical Framework

The proposed numerical framework is predicated on a fully coupled pressure-based algorithm with a conservative finite-volume discretization of the governing equations and a collocated variable arrangement [18], with a spatial and temporal convergence that is consistent with the employed second-order schemes [18,30]. The momentum-weighted interpolation used to define the fluxes at cell faces is the only source of numerical dissipation, leading to a small error in kinetic energy that converges with third order under mesh refinement [18,35]. The numerical framework presented in the following assumes a rectilinear mesh, without skewness or non-orthogonality. Corrections for mesh skewness and non-orthogonality may be included in the presented discretization following the work of Denner et al. [18].

4.1. Finite-Volume Discretization

The BDF2 scheme, also known as Second-Order Backward Euler scheme, is used to discretize the transient terms of the governing equations, which is defined for the transient term of a general fluid variable, $\varphi$, at cell P as [18]

$${\displaystyle {\int \int \int }_V}{\displaystyle \frac{\partial \varphi }{\partial t}}\mathrm{d}V\approx {\left.\frac{\partial \varphi }{\partial t}\right|}_P{V}_P=\left[\left({\displaystyle \frac{1}{\Delta t_1}}+{\displaystyle \frac{1}{\Delta \tau }}\right){\varphi }_P-\left({\displaystyle \frac{1}{\Delta t_1}}+\frac{1}{\Delta t_2}\right){\varphi }_P^{(t-\Delta t_1)}+{\displaystyle \frac{\Delta t_1}{\Delta t_2\Delta \tau }}{\varphi }_P^{(t-\Delta \tau )}\right]V_P,$$

(10)

with $\Delta \tau =\Delta t_1+\Delta t_2$, where V is the volume of the mesh cell, $\Delta t_1$ is the current time-step and $\Delta t_2$ is the previous time-step. The superscripts $(t-\Delta t_1)$ and $(t-\Delta \tau )$ denote the values of the previous and previous-previous time-level, respectively. For consistency, the transient terms of both the continuity Equation (1) and the momentum Equation (2) are discretized with the same scheme.

Applying the divergence theorem and the midpoint rule, the advection terms of the governing equations are discretized as

$${\displaystyle {\int \int \int }_V}\frac{\partial \rho u_i\varphi }{\partial x_i}\mathrm{d}V\approx \sum _f{\tilde{\rho }}_f{\vartheta }_f{\tilde{\varphi }}_f{A}_f,$$

(11)

where subscript f denotes the mesh face shared by cell P and its neighbor cell Q, as illustrated in Figure 1, ${\vartheta }_f={\mathit{u}}_f·{\mathit{n}}_f$ is the advecting velocity at face f, presented in detail in Section 4.2, $A_f$ is the area of mesh face f and ${\mathit{n}}_f$ is the unit normal vector of face f, which points out of cell P. The advected variable ${\tilde{\varphi }}_f$ is interpolated using a TVD formulation, given as [36]

$${\tilde{\varphi }}_f\approx {\varphi }_U+{\xi }_f\frac{|x_f-x_U|}{\Delta x_f}\left({\varphi }_D-{\varphi }_U\right),$$

(12)

where ${\varphi }_U$ and ${\varphi }_D$ are the values of the advected variable at the upwind and downwind cells, respectively, $x_U$ is the position of the center of the upwind cell U, $x_f$ is the position of the center of face f, $\Delta x_f$ is the distance between the cell centers adjacent to face f, and ${\xi }_f$ is the flux limiter.

Applying the divergence theorem, the viscous stress term of the momentum Equation (2) is given as

$$\begin{array}{cc}\hfill {\displaystyle {\int \int \int }_V}\frac{\partial {\tau }_{ji}}{\partial x_i}\mathrm{d}V& \approx \sum _f{\mu }_f^{*}\left({\left.\frac{\partial u_j}{\partial x_i}\right|}_f+{\left.\overline{\frac{\partial u_i}{\partial x_j}}\right|}_f-\frac{2}{3}{\left.\overline{\frac{\partial u_k}{\partial x_k}}\right|}_f\right)n_{i,f}{A}_f,\hfill \end{array}$$

(13)

where $\overline{\square }$ denotes linear interpolation from the adjacent cell centers and ${\mu }_f^{*}$ is the face viscosity determined by harmonic interpolation from the adjacent cell centers. The face-centered gradient of the primary velocity component, $u_j$, is discretized as

$${\left.\frac{\partial u_j}{\partial x_i}\right|}_f{n}_{i,f}\approx \frac{u_{j,Q}-u_{j,P}}{\Delta x_f}.$$

(14)

4.2. Advecting Velocity

The advecting velocity, ${\vartheta }_f={\mathit{u}}_f·{\mathit{n}}_f$, is discretized based on a momentum-weighted interpolation (MWI) [37]. The MWI provides a robust coupling between pressure and velocity for flows at low Mach numbers [35] and avoids pressure-velocity decoupling as a result of the applied collocated variable arrangement. The advecting velocity ${\vartheta }_f$ at face f is given as [35]

$$\begin{array}{cc}\hfill {\vartheta }_f={\overline{\mathit{u}}}_f·{\mathit{n}}_f& -{\widehat{d}}_f\left[\frac{p_Q-p_P}{\Delta x_f}-{\rho }_f^{*}\left(\frac{1-l_{Pf}}{{\rho }_P}{\left.\frac{\partial p}{\partial x_i}\right|}_P+\frac{l_{Pf}}{{\rho }_Q}{\left.\frac{\partial p}{\partial x_i}\right|}_Q\right)n_{i,f}\right]\hfill \\ & +{\widehat{d}}_f\frac{{\rho }_f^{*(t-\Delta t_1)}}{\Delta t_1}\left({\vartheta }_f^{(t-\Delta t_1)}-{\overline{\mathit{u}}}_{i,f}^{(t-\Delta t_1)}·{\mathit{n}}_{i,f}\right),\hfill \end{array}$$

(15)

where the face density ${\rho }_f^{*}$ is determined by harmonic interpolation from the cell centers adjacent to face f. The coefficient ${\widehat{d}}_f$ is defined based on the applied time-step as well as the velocity coefficients associated with the advection and shear stress terms of the discretized momentum equations [35]. The density weighting applied to the cell-centered pressure gradient has been shown to yield robust results for flows with large and abrupt changes in density [35], demonstrated for incompressible interfacial flows with a density ratio of up to ${10}^{24}$ [26,38]. The transient term of Equation (15) is critical for a correct representation of pressure waves and ensures that the MWI is time-step independent [35].

4.3. Discretized Governing Equations

Applying the discretization schemes described above, the discretized continuity Equation (1) is given at cell P as

$${\left.{\displaystyle \frac{\partial \rho }{\partial t}}\right|}_P{V}_P+\sum _f{\tilde{\rho }}_f^{(n+1)}{\vartheta }_f^{(n+1)}{A}_f=0$$

(16)

and the discretized momentum Equation (2) are given at cell P as

$$\begin{array}{cc}\hfill {\left.{\displaystyle \frac{\partial \rho u_j}{\partial t}}\right|}_P{V}_P& +\sum _f{\tilde{\rho }}_f^{(n+1)}{\vartheta }_f^{(n+1)}{\tilde{u}}_{j,f}^{(n+1)}{A}_f=-\sum _f{\overline{p}}_f^{(n+1)}{n}_{j,f}{A}_f\hfill \\ & +\sum _f{\mu }_f^{*}\left(\frac{u_{j,Q}^{(n+1)}-u_{j,P}^{(n+1)}}{\Delta x_f}+{\left.\overline{\frac{\partial u_i}{\partial x_j}}\right|}_f^{\left(n\right)}{n}_{i,f}-\frac{2}{3}{\left.\overline{\frac{\partial u_k}{\partial x_k}}\right|}_f^{\left(n\right)}{n}_{i,f}\right)A_f,\hfill \end{array}$$

(17)

where superscript n is the iteration counter of the nonlinear iterations conducted to solve the discretized governing equations in each time-step, see Section 4.4. To this end, superscript $(n+1)$ denotes the solution that is sought implicitly and superscript $\left(n\right)$ stands for the most recent available solution. A Newton linearization is applied to linearize both governing equations [29]. The linearization of the advection term of the continuity equation as well as the transient term of the momentum equations is defined as

$${\varphi }_1^{(n+1)}{\varphi }_2^{(n+1)}\approx {\varphi }_1^{\left(n\right)}{\varphi }_2^{(n+1)}+{\varphi }_1^{(n+1)}{\varphi }_2^{\left(n\right)}-{\varphi }_1^{\left(n\right)}{\varphi }_2^{\left(n\right)}$$

(18)

and the linearization of the advection term of the momentum equations is defined as

$${\varphi }_1^{(n+1)}{\varphi }_2^{(n+1)}{\varphi }_3^{(n+1)}\approx {\varphi }_1^{\left(n\right)}{\varphi }_2^{\left(n\right)}{\varphi }_3^{(n+1)}+{\varphi }_1^{\left(n\right)}{\varphi }_2^{(n+1)}{\varphi }_3^{\left(n\right)}+{\varphi }_1^{(n+1)}{\varphi }_2^{\left(n\right)}{\varphi }_3^{\left(n\right)}-2{\varphi }_1^{\left(n\right)}{\varphi }_2^{\left(n\right)}{\varphi }_3^{\left(n\right)},$$

(19)

where ${\varphi }_1$, ${\varphi }_2$ and ${\varphi }_3$ are generic fluid variables. This linearization facilitates an implicit treatment of all terms in the governing equations that depend on pressure and velocity, which has previously been shown to improve the stability and performance of the solution algorithm for flows in all Mach number regimes [29,39,40]. It also provides the implicit pressure-velocity coupling necessary to robustly predict flows with low Mach numbers.

The implicit treatment of the density at cell centers is given by reformulating Equation (8) as an implicit function of the cell-centered pressure, given as

$${\rho }_P^{(n+1)}=\frac{K{\left(p_P^{\left(n\right)}+\mathsf{\Pi }\right)}^{\mathsf{\Gamma }-1}}{1+bK{\left(p_P^{\left(n\right)}+\mathsf{\Pi }\right)}^{\mathsf{\Gamma }}}\left(p_P^{(n+1)}+\mathsf{\Pi }\right)$$

(20)

and the implicit treatment of the advecting velocity is given by the semi-implicit formulation

$$\begin{array}{cc}\hfill {\vartheta }_f={\overline{\mathit{u}}}_f^{(n+1)}·{\mathit{n}}_f& -{\widehat{d}}_f\left[\frac{p_Q^{(n+1)}-p_P^{(n+1)}}{\Delta x_f}-{\rho }_f^{*\left(n\right)}\left(\frac{1-l_{Pf}}{{\rho }_P^{\left(n\right)}}{\left.\frac{\partial p}{\partial x_i}\right|}_P^{\left(n\right)}+\frac{l_{Pf}}{{\rho }_Q^{\left(n\right)}}{\left.\frac{\partial p}{\partial x_i}\right|}_Q^{\left(n\right)}\right)n_{i,f}\right]\hfill \\ & +{\widehat{d}}_f\frac{{\rho }_f^{*(t-\Delta t_1)}}{\Delta t_1}\left({\vartheta }_f^{(t-\Delta t_1)}-{\overline{\mathit{u}}}_{i,f}^{(t-\Delta t_1)}·{\mathit{n}}_{i,f}\right).\hfill \end{array}$$

(21)

4.4. Solution Procedure

The linear system of discretized governing equations, $\mathcal{A}\mathit{\varphi }=\mathit{b}$, is solved simultaneously for the pressure p and the velocity vector $\mathit{u}\equiv {(u,v,w)}^T$. Considering a three-dimensional computational mesh with N cells, this linear system of equations then follows as

$$\left(\begin{array}{cccc}{\mathcal{A}}^{\rho ,p}& {\mathcal{A}}^{\rho ,u}& {\mathcal{A}}^{\rho ,v}& {\mathcal{A}}^{\rho ,w}\\ {\mathcal{A}}^{\rho u,p}& {\mathcal{A}}^{\rho u,u}& {\mathcal{A}}^{\rho u,v}& {\mathcal{A}}^{\rho u,w}\\ {\mathcal{A}}^{\rho v,p}& {\mathcal{A}}^{\rho v,u}& {\mathcal{A}}^{\rho v,v}& {\mathcal{A}}^{\rho v,w}\\ {\mathcal{A}}^{\rho w,p}& {\mathcal{A}}^{\rho w,u}& {\mathcal{A}}^{\rho w,v}& {\mathcal{A}}^{\rho w,w}\end{array}\right)·\left(\begin{array}{c}{\mathit{\varphi }}^p\\ {\mathit{\varphi }}^u\\ {\mathit{\varphi }}^v\\ {\mathit{\varphi }}^w\end{array}\right)=\left(\begin{array}{c}{\mathit{b}}^{\rho }\\ {\mathit{b}}^{\rho u}\\ {\mathit{b}}^{\rho v}\\ {\mathit{b}}^{\rho w}\end{array}\right),$$

(22)

where ${\mathcal{A}}^{\eta ,\chi }$ are the coefficient submatrices of size $N\times N$, where $\eta$ denotes the conserved quantity and $\chi$ denotes the primary solution variable associated with a given governing equation. The subvectors ${\mathit{\varphi }}^{\chi }$ with length N contain the solution of the primary solution variable $\chi$ and the subvectors ${\mathit{b}}^{\eta }$ with length N contain all contributions from boundary conditions, previous nonlinear iterations as well as previous time-levels. The solution procedure conducts nonlinear iterations to solve the system of discretized Equations (22), in which (22) is preconditioned by the Block-Jacobi preconditioner and then solved using the BiCGSTAB solver of the software library PETSc [41,42,43], as discussed in detail by Denner et al. [18].

5. Interface Treatment

The discretized governing equations are extended to interfacial flows using an interface advection method (see Section 5.1) in conjunction with an appropriate definition of the fluid properties (see Section 5.2). In order to represent the two interacting fluids discretely, the indicator function $\zeta$ translates into a color function $\psi$, defined for cell P as

$${\psi }_P=\frac{1}{V_P}{\displaystyle {\int \int \int }_V}\zeta \mathrm{d}V.$$

(23)

5.1. Interface Advection

To advect the fluid interface between two fluids, Equation (5) is integrated in each cell, which gives the following transport equation for the color function $\psi$, Equation (23),

$$\frac{\partial \psi }{\partial t}+\frac{\partial u_i\psi }{\partial x_i}-\psi \frac{\partial u_i}{\partial x_i}=0,$$

(24)

which is discretized using an algebraic VOF method [30]. Discretizing the transient term using the Crank–Nicolson scheme, the semi-discretized form of Equation (24) is given as

$$\frac{{\psi }_P-{\psi }_P^{(t-\Delta t_{\psi })}}{\Delta t_{\psi }}{V}_P+\sum _f\frac{{\psi }_f+{\psi }_f^{(t-\Delta t_{\psi })}}{2}{\vartheta }_f{A}_f-\frac{{\psi }_P+{\psi }_P^{(t-\Delta t_{\psi })}}{2}\sum _f{\vartheta }_f{A}_f=0,$$

(25)

where $\Delta t_{\psi }$ is the applied VOF time-step. The same advecting velocity ${\vartheta }_f$ that is used for all advection terms of the discretized governing equations is applied to advect the color function. The face value ${\psi }_f$ is discretized using the CICSAM scheme [44], which accounts for the available flux volume as well as the orientation of the interface. To retain a sharp interface, a relatively small time-step $\Delta t_{\psi }$ is used, satisfying a Courant number of ${\mathrm{Co}}_{\psi }=\Delta t_{\psi }\left|\mathit{u}\right|/\Delta x\le 0.05$, and a dual time-stepping method [45] is applied for a better computational performance. The time-step $\Delta t_{\psi }$ applied to solve the VOF advection equation is generally smaller than the time-step $\Delta t_1$ applied to solve Equation (22). The fluid time-step $\Delta t_1$, thus, has to be an integer multiple of the VOF time-step $\Delta t_{\psi }$ to ensure a consistent advection of the color function.

This algebraic VOF method captures volume changes of the bulk phases with second-order accuracy [30] and has been used successfully for studies of bubble dynamics in compressible [15,30] and incompressible [46] interfacial flows. Nevertheless, the presented pressure-based algorithm is not limited to the employed algebraic VOF method; other methods to represent the bulk phases or advect the interface may equally be applied.

5.2. Fluid Properties

The treatment of density at fluid interfaces is based on the acoustically conservative interface discretization (ACID) method [30]. The primary assumption of the ACID method is that all cells in the finite-volume stencil of a given cell are assigned the same color function value, i.e., the color function is assumed to be constant in the entire finite-volume stencil. Density, which is discontinuous at the interface, is then evaluated based on this locally constant color function field, which recovers the contact discontinuity represented by the interface [33,47]. Denner et al. [30] reported robust and accurate results for acoustic and shock waves in interfacial flows, supporting the notion that the interface discretization indeed conserves the acoustic features of the flow.

Under the assumptions made by the ACID method, the density interpolated to face f of cell P is given as

$${\tilde{\rho }}_f={\rho }_U^{🟉}+{\xi }_f\frac{|x_f-x_U|}{\Delta x_f}\left({\rho }_D^{🟉}-{\rho }_U^{🟉}\right).$$

(26)

The density ${\rho }_U$ at the upwind cell U of face f and the density ${\rho }_D$ at the downwind cell D of face f are defined based on the color function value of cell P, ${\psi }_P$, as

$$\begin{array}{cc}\hfill {\rho }_U^{🟉}& =(1-{\psi }_P){\rho }_{\mathrm{a},U}+{\psi }_P{\rho }_{\mathrm{b},U}\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\rho }_D^{🟉}& =(1-{\psi }_P){\rho }_{\mathrm{a},D}+{\psi }_P{\rho }_{\mathrm{b},D},\hfill \end{array}$$

(28)

where the partial densities ${\rho }_{\mathrm{a}}$ and ${\rho }_{\mathrm{b}}$ associated with fluids $\mathrm{a}$ and $\mathrm{b}$ are given by Equation (8). The densities at previous time-levels are evaluated similarly, based on ${\psi }_P$, following as

$${\rho }_P^{(t-\Delta t_1)}=(1-{\psi }_P){\rho }_{\mathrm{a},P}^{(t-\Delta t_1)}+{\psi }_P{\rho }_{\mathrm{b},P}^{(t-\Delta t_1)}$$

(29)

and analogously for ${\rho }_P^{(t-\Delta \tau )}$.

In order to treat the density implicitly, as given in Equation (20), the coefficient of the implicit pressure value is defined as

$$\frac{K{\left(p^{\left(n\right)}+\mathsf{\Pi }\right)}^{\mathsf{\Gamma }-1}}{1+bK{\left(p^{\left(n\right)}+\mathsf{\Pi }\right)}^{\mathsf{\Gamma }}}=(1-\psi )\frac{K_{\mathrm{a}}{\left(p^{\left(n\right)}+{\mathsf{\Pi }}_{\mathrm{a}}\right)}^{{\mathsf{\Gamma }}_{\mathrm{a}}-1}}{1+b_{\mathrm{a}}{K}_{\mathrm{a}}{\left(p^{\left(n\right)}+{\mathsf{\Pi }}_{\mathrm{a}}\right)}^{{\mathsf{\Gamma }}_{\mathrm{a}}}}+\psi \frac{K_{\mathrm{b}}{\left(p^{\left(n\right)}+{\mathsf{\Pi }}_{\mathrm{b}}\right)}^{{\mathsf{\Gamma }}_{\mathrm{b}}-1}}{1+b_{\mathrm{b}}{K}_{\mathrm{b}}{\left(p^{\left(n\right)}+{\mathsf{\Pi }}_{\mathrm{b}}\right)}^{{\mathsf{\Gamma }}_{\mathrm{b}}}}$$

(30)

and the pressure constant $\mathsf{\Pi }$ as

$$\frac{K{\left(p^{\left(n\right)}+\mathsf{\Pi }\right)}^{\mathsf{\Gamma }-1}}{1+bK{\left(p^{\left(n\right)}+\mathsf{\Pi }\right)}^{\mathsf{\Gamma }}}\mathsf{\Pi }=(1-\psi )\frac{K_{\mathrm{a}}{\left(p^{\left(n\right)}+{\mathsf{\Pi }}_{\mathrm{a}}\right)}^{{\mathsf{\Gamma }}_{\mathrm{a}}-1}}{1+b_{\mathrm{a}}{K}_{\mathrm{a}}{\left(p^{\left(n\right)}+{\mathsf{\Pi }}_{\mathrm{a}}\right)}^{{\mathsf{\Gamma }}_{\mathrm{a}}}}{\mathsf{\Pi }}_{\mathrm{a}}+\psi \frac{K_{\mathrm{b}}{\left(p^{\left(n\right)}+{\mathsf{\Pi }}_{\mathrm{b}}\right)}^{{\mathsf{\Gamma }}_{\mathrm{b}}-1}}{1+b_{\mathrm{b}}{K}_{\mathrm{b}}{\left(p^{\left(n\right)}+{\mathsf{\Pi }}_{\mathrm{b}}\right)}^{{\mathsf{\Gamma }}_{\mathrm{b}}}}{\mathsf{\Pi }}_{\mathrm{b}},$$

(31)

with the polytropic constants $K_{\mathrm{a}}$ and $K_{\mathrm{b}}$ given by Equation (9) based on the reference pressure and reference density of the respective bulk phase.

The viscosity $\mu$ in cell P is defined by a linear interpolation based on the color function, given as

$$\begin{array}{c}\hfill \mu =(1-\psi ){\mu }_{\mathrm{a}}+\psi {\mu }_{\mathrm{b}}.\end{array}$$

(32)

6. Results

The proposed algorithm is validated for the simulation of acoustic cavitation using test-cases representative of this kind of flow. The propagation of acoustic waves presented in Section 6.1 tests the accurate treatment of acoustic effects and the Rayleigh collapse in Section 6.2 tests the correct prediction of the entire pressure-driven bubble dynamics. Lastly, the wall-bounded cavitation of a bubble for a broad range of bubble stand-off distances in Section 6.3 scrutinizes the reliable prediction of complex cavitation events, also in comparison to experiments. For this validation, four different fluids are considered, with the fluid properties given in

Table 1. Properties of the considered fluids.

Fluid	$\mathit{\kappa }$	b   [m${}^3$  kg${}^{-1}$]	$\mathsf{\Pi }$ [Pa]
Air	$1.400$	0	0
JA2 propellant gas [48]	$1.225$	$1.00\times {10}^{-3}$	0
Water Tait [23]	$7.150$	0	$3.046\times {10}^8$
Water NASG [32]	$1.187$	$6.61\times {10}^{-4}$	$7.028\times {10}^8$

. While the Tait model of water is frequently used to simulate cavitation events in water, using the full NASG model of water provides an alternative description. The JA2 propellant gas is considered to be an alternative gas, to demonstrate the prediction of acoustic effects in a Noble–Abel gas.

6.1. Acoustic Waves

Reproducing the interaction between acoustic waves and fluid interfaces correctly is a basic requirement for the accurate prediction of acoustic cavitation. Given a small perturbation to the flow, with the amplitude of the velocity perturbation $\Delta u_0\ll a_0$, the resulting acoustic wave propagates with the speed of sound $a_0$, given for a general NASG fluid as

$$a_0=\sqrt{\kappa \frac{p_0+\mathsf{\Pi }}{{\rho }_0(1-b{\rho }_0)}},$$

(33)

with a pressure wave of amplitude $\Delta p_0={\rho }_0{a}_0\Delta u_0$ [47]. For all considered cases, the unperturbed flow velocity is $u_0=1\mathrm{m}{\mathrm{s}}^{-1}$ and the ambient pressure is $p_0={10}^5\mathrm{Pa}$. The computational domain is represented with a mesh spacing of $\Delta x=2\times {10}^{-3}\mathrm{m}$.

First, the propagation of acoustic waves in a single-phase flow is considered. The acoustic waves are generated by a periodically changing velocity at the domain inlet, which is defined as $u_{\mathrm{in}}=u_0+\Delta u_{0}sin\left(2\pi ft\right)$, where $\Delta u_0=0.01u_0$ is the amplitude of the velocity perturbation and f is the excitation frequency given in

Table 2. The excitation frequency f of the acoustic waves, the reference density ${\rho }_0$ and the associated speed of sound $a_0$ at the reference pressure of $p_0={10}^5\mathrm{Pa}$, the wavelength ${\lambda }_0$ and pressure amplitude $\Delta p_0$ of the acoustic waves given by linear acoustic theory, and the wavelength $\lambda$ and  the pressure amplitude $\Delta p$ of the acoustic waves predicted by the proposed algorithm.

Fluid	f  [s${}^{-1}$]	${\mathit{\rho }}_0$  [kg  m${}^{-3}$]	${\mathit{a}}_0$  [m  s${}^{-1}$]	${\mathit{\lambda }}_0$ [m]	$\mathbf{\Delta }{\mathit{p}}_0$ [Pa]	$\mathit{\lambda }$ [m]	$\mathbf{\Delta }\mathit{p}$ [Pa]
Air	1750	$1.157$	$347.85$	$0.199$	$4.025$	$0.199$	$4.019$
JA2 propellant gas [48]	1750	$0.997$	$350.70$	$0.200$	$3.496$	$0.201$	$3.492$
Water Tait [23]	7500	1000	$1476.0$	$0.197$	$14,760$	$0.198$	$14,741$
Water NASG [32]	7500	1000	$1568.8$	$0.209$	$15,688$	$0.210$	$15,673$

. Figure 2 shows the pressure profiles of acoustic waves in each fluid shortly before the waves reach the domain outlet at $x=1.0\mathrm{m}$. The wavelength $\lambda$ and pressure amplitude $\Delta p$ predicted by the proposed algorithm for these acoustic waves are given in Table 2. For all cases $\lambda$ and $\Delta p$ are in excellent agreement with the theoretical values ${\lambda }_0$ and $\Delta p_0$ according to linear acoustic theory, also given in Table 2.

The propagation of a single acoustic wave in gas-liquid flows is simulated to demonstrate the accurate prediction of the interaction of acoustic waves with fluid interfaces. The acoustic wave is initiated in the left phase by a Gaussian pressure pulse as

$$p\left(x\right)=p_0+\Delta p_0{e}^{-\frac{1}{2}{\left(\frac{x-x_0}{\sigma }\right)}^2},$$

(34)

where $\Delta p_0$ and $x_0$ are the initial amplitude and the initial position of the pressure pulse, respectively. Based on linear acoustic theory [47], the pressure pulse reflected in the left phase at the fluid interface should have a pressure amplitude of

$$\Delta p_{\mathrm{L},0}^{\mathrm{reflected}}={\displaystyle \frac{\Delta p_{\mathrm{L},0}^{\mathrm{incident}}}{{\displaystyle \frac{2Z_{\mathrm{R}}}{Z_{\mathrm{R}}-Z_{\mathrm{L}}}}-1}},$$

(35)

with the pressure amplitude of the incident pulse given as $\Delta p_{\mathrm{L},0}^{\mathrm{incident}}=\Delta p_0$, and the pressure amplitude of the pulse transmitted to the right phase should be $\Delta p_{\mathrm{R},0}^{\mathrm{transmitted}}=p_{\mathrm{L},0}^{\mathrm{incident}}+\Delta p_{\mathrm{L},0}^{\mathrm{reflected}}$, where subscripts L and R denote the left and right phases, respectively, and $Z=\rho a$ is the characteristic acoustic impedance. The pressure profiles of the pulses are shown in Figure 3 for two air–water systems, after the pressure pulse interacted with the fluid interface. Figure 3a shows the pressure pulse, initialized in the gas phase with $\Delta p_0=10\mathrm{Pa}$ and $\sigma =0.03$ at $x_0=0.2\mathrm{m}$, in an air–water system, with water described using the properties of Water NASG given in Table 1. The pressure amplitude of the reflected and transmitted pulses are both in excellent agreement with linear acoustic theory. Similarly, the pressure amplitudes of the reflected and transmitted pulses are also in excellent agreement with linear acoustic theory for the case shown in Figure 3b, where the pressure pulse is initialized in the liquid phase with $\Delta p_0=1000\mathrm{Pa}$ and $\sigma =0.1$ at $x_0=0.6\mathrm{m}$, and with water described using the properties of Water Tait given in Table 1.

6.2. Rayleigh Collapse

The behavior of a spherical bubble collapsing as a result of an initial overpressure in the liquid compared to the pressure in the gas bubble is considered, the so-called Rayleigh collapse [49], to scrutinize the prediction of pressure-driven bubble dynamics by the proposed algorithm. The Gilmore model [9], which is an extension of the classical RP model, for the dynamic behavior of a spherical gas bubble in a compressible liquid is used as a reference solution. The solution of the Gilmore model is, therefore, regarded as the exact solution of the dynamic bubble behavior under consideration. The dissipation due to both viscous stresses and acoustic radiation is included in these simulations, leading to a complex interplay of the hydrodynamics, e.g., viscous dissipation, and thermodynamics, e.g., energy dissipation due to the emission of the shock wave formed at the end of the collapse, which directly influences the evolution of the bubble radius. In particular, the accurate prediction of the bubble radius after multiple collapse-expansion cycles using finite-volume methods is known to be difficult to achieve [23,50].

The initial pressure inside the bubble is $p_{\mathrm{b}}=4\times {10}^3\mathrm{Pa}$ and the initial pressure of the liquid is $p_{\mathrm{l}}\left(R\right)=p_{\infty }+(p_{\mathrm{b}}-p_{\infty })R_0/R$, with $p_{\infty }={10}^5\mathrm{Pa}$. The gas is modeled as air by the ideal-gas model with the properties given in Table 1 and a reference density of ${\rho }_{0,\mathrm{air}}=1.2\mathrm{kg}{\mathrm{m}}^{-3}$, and a viscosity of ${\mu }_{\mathrm{air}}=1.82\times {10}^{-5}\mathrm{Pa}\mathrm{s}$. The liquid is taken to be water described by the Tait model with the properties given for Water Tait in Table 1 and a reference density ${\rho }_{0,\mathrm{water}}=1000\mathrm{kg}{\mathrm{m}}^{-3}$, and a viscosity of ${\mu }_{\mathrm{water}}={10}^{-3}\mathrm{Pa}\mathrm{s}$. The reference pressure is $p_0=p_{\infty }$ for both gas and liquid. Figure 4 shows the dimensionless bubble radius $R/R_0$, normalized by the initial radius $R_0$ (which also corresponds to the maximum radius in this case), as a function of the dimensionless time $t/t_{\mathrm{c}}$, normalized by the Rayleigh collapse time $t_{\mathrm{c}}=0.915R_0\sqrt{{\rho }_{0,\mathrm{water}}/p_{\infty }}$, obtained with different spatial and temporal resolutions, compared against the solution obtained using the Gilmore model. Excellent agreement between the evolution of the bubble radius R predicted by the proposed algorithm and the Gilmore solution is observed in Figure 4, even after multiple collapse-expansion cycles, with the simulation results converging to the solution of the Gilmore model for sufficiently large spatial and temporal resolution.

6.3. Wall-Bounded Cavitation

The wall-bounded collapse of a laser-induced cavitation bubble is simulated to test the fidelity of the proposed algorithm in predicting a complex pressure-driven cavitation process. The bubble is initialized with a radius of $R_0=50\mathsf{\mu }\mathrm{m}$ and a pressure $p_{\mathrm{b}}>p_{\infty }$ to achieve a maximum radius of the upper hemisphere of the bubble of $R_{\max }\approx 390\mathsf{\mu }\mathrm{m}$. The bubble is situated in water with an ambient pressure of $p_{\infty }={10}^5\mathrm{Pa}$ at a distance ${\ell }_0$ from the wall, as illustrated in Figure 5a. The reference density at $p_0=p_{\infty }$ for air and water is ${\rho }_{0,\mathrm{air}}=1.2\mathrm{kg}{\mathrm{m}}^{-3}$ and ${\rho }_{0,\mathrm{water}}=1000\mathrm{kg}{\mathrm{m}}^{-3}$, respectively. Water is modeled using the properties of Water NASG given in Table 1 and has a viscosity of ${\mu }_{\mathrm{water}}={10}^{-3}\mathrm{Pa}\mathrm{s}$. Air is modeled as an ideal gas, using the properties given in Table 1, and has a viscosity of ${\mu }_{\mathrm{air}}=1.82\times {10}^{-5}\mathrm{Pa}\mathrm{s}$. Apart from the bottom wall illustrated explicitly in Figure 5, all boundaries of the computational domain are situated at a distance of $0.25\mathrm{m}$ from the bubble inception site, so that acoustic waves cannot reach the boundaries within the time span considered in the simulations. The computational domain is resolved with a mesh spacing of $\Delta x=2\mathsf{\mu }\mathrm{m}$ and a gradually refined mesh resolution near the wall, with a minimum mesh spacing of $\Delta x_{\min }=0.05\mathsf{\mu }\mathrm{m}$. The mesh is rapidly coarsened outside the region that the bubble occupies. The time-step is chosen to satisfy the Courant number $\mathrm{Co}=\Delta t\left|\mathbf{u}\right|/\Delta x\le 0.8$.

Figure 5b–d show the instantaneous bubble shape and contours of the velocity component perpendicular to the wall of a bubble with $\gamma ={\ell }_0/R_{\max }=0.59$ and $R_{\max }=388\mathsf{\mu }\mathrm{m}$ at times $t=39.0\mathsf{\mu }\mathrm{s}$, $t=70.5\mathsf{\mu }\mathrm{s}$ and $t=85.0\mathsf{\mu }\mathrm{s}$, respectively. As the bubble initially grows and then collapses, a thin liquid film forms between the bubble and the wall, which can be clearly seen in Figure 5b–d and which has also previously been observed and measured in experiments [51]. Furthermore, the characteristic fast-moving liquid jet directed towards the wall that pierces the bubble during its collapse can be seen in Figure 5d. Figure 6a shows the minimum thickness, ${\ell }_{\min }$, of the liquid film separating the bubble and the wall at the time when the jet pierces through the bubble, i.e., the thickness of the small gap between the bubble and the wall seen in Figure 5d, predicted by the simulations for different dimensionless initial stand-off distances $\gamma$, in comparison with the experimental measurements of Reuter and Kaiser [51]. The fit of the correlation has a coefficient of determination of $R^2=0.94$ based on 91 experiments for bubbles with $R_{\max }=385-410\mathsf{\mu }\mathrm{m}$ in the range $\gamma =0.47-1.07$ [51]. Excellent agreement between the simulations and the experiments is observed, demonstrating the accuracy with which the proposed algorithm can predict such a complex cavitation process. The evolution of the peak shear stress and peak pressure at the wall are shown in Figure 6b,c, respectively, for selected initial stand-off distances. The pressure signal exhibits two peaks during the collapse of the bubble [52]; the first peak is caused by the pressure pulse emitted as the liquid jet impinges on the lower side of the bubble interface and the second peak is associated with the shock wave emitted at the end of the collapse, when the bubble (at this point having a toroidal shape) assumes its minimum volume. The results are in good agreement with previous studies regarding the order of magnitude of the generated wall shear stress [11,12] and the pressure amplitude of the shock wave formed as a result of the collapse [53,54] of a laser-induced bubble, with wall shear stresses exceeding $100\mathrm{kPa}$ for bubbles collapsing close to a wall.

7. Conclusions

A fully coupled pressure-based algorithm for the simulation of the acoustic cavitation of bubbles in polytropic gas-liquid systems has been proposed. The algorithm is based on a conservative finite-volume discretization with a collocated variable arrangement [18], whereby the discretized governing conservation laws are solved simultaneously in a single linear system of equations for pressure and velocity, with density described by a polytropic assumption using the Noble-Abel stiffened-gas model [32]. The interface separating the interacting bulk phases is captured by a state-of-the-art algebraic VOF method [30].

The proposed numerical framework has been validated by representative test-cases, including the propagation and interface interaction of acoustic waves, as well as the Rayleigh collapse and the wall-bounded cavitation of a bubble. The propagation of acoustic waves demonstrates the accurate prediction of acoustic effects, whereas the excellent comparison of the bubble radius of the Rayleigh collapse against the Gilmore model demonstrates the correct prediction of the entire pressure-driven bubble dynamics. In particular, the demonstrated accurate prediction of the bubble radius after multiple collapse-expansion cycles compared to the theoretical solution, in this case given by the Gilmore model, is known to be difficult to achieve [23,50]. The considered wall-bounded cavitation of a bubble shows the reliable prediction of complex cavitation events. Especially the predicted thickness of the thin liquid film, which separates the bubble from the wall, is in very good agreement with recent experimental measurements [51] for a broad range of bubble stand-off distances, further establishing the reliable prediction of acoustic cavitation by the proposed algorithm.