# An Assessment of CFD Cavitation Models Using Bubble Growth Theory and Bubble Transport Modeling

Michael P Kinzel, Robert F Kunz, Jules W Lindau

► To cite this version:

Michael P Kinzel, Robert F Kunz, Jules W Lindau. An Assessment of CFD Cavitation Models Using Bubble Growth Theory and Bubble Transport Modeling. 17th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery (ISROMAC2017), Dec 2017, Maui, United States. hal-02419952

HAL Id: hal-02419952

https://hal.science/hal-02419952

Submitted on 19 Dec 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An Assessment of CFD Cavitation Models Using Bubble Growth Theory and Bubble Transport Modeling

Michael P. Kinzel<fnref n="1" />*, Robert F. Kunz<fnref n="2" />, Jules W. Lindau <fnref n="1" />

<figure><img src="image_2.png" /></figure>

ISROMAC 2017

International Symposium on Transport Phenomena and Dynamics of Rotating Machinery

Maui, Hawaii

December 16-21, 2017

## Abstract

This effort investigates two approaches to cavitation modeling that are relevant to computational fluid dynamics (CFD). The two approaches include (1) reformulating the cavitation models and (2) exploring the impact of liquid-vapor slip. The first aspect of the paper revisits cavitation model formulations with respect to the Rayleigh-Plesset Equation (RPE). The approach reformulates the cavitation model using analytic solutions to the RPE. The benefit of this reformulation is displayed through maintaining model sensitivities similar to RPE, whereas the standard models fail these tests. In addition, the model approach is extended beyond standard homogenous models to a two-fluid model framework that includes bubble slip. The results indicate a significant impact of slip on the predicted cavitation solution, suggesting that inclusion of such modeling can potentially improve CFD cavitation models. Overall, the results of this effort point to various aspects that should be considered in future CFD-modeling efforts that aim to model cavitation.

## Keywords

CFD – Cavitation Modeling – Approximate Bubble Dynamics – Bubble Transport

## INTRODUCTION

Cavitation relates to marine applications, pump designs, rockets, and a variety of other liquid flows. In general, this is an undesired physical process that is highly non-linear and strongly impacts performance. For these reasons, modeling cavitation using computational fluid dynamics (CFD) has become important for the design of complex systems. Because of this, maturing physically accurate models for cavitation remains useful.

In the context of CFD, there are various approaches to cavitation modeling. A common approach is based on a volume of fluids (VOF)-like approach, coupled with a model for the exchange from a liquid to vapor phase (evaporation) and vice versa (condensation). These models are commonly referred to as homogenous multiphase models [1, 2, 3, 4]. These approaches have all improved significantly in their reliability over the years, but still pose challenges in complexity, mesh requirements, and difficult-to-ascertain empirical factors. Such a model aims to directly conserve and transport vapor mass and uses source terms to model the evaporation and condensation processes. A vapor-mass conservation equation used is given by [1, 2]:

$$ \frac{\partial(\rho_v \alpha_v)}{\partial t} + \frac{\partial(\rho_v \alpha_v u_i)}{\partial x_i} = \dot{m}_{evap} - \dot{m}_{cond} \quad (1) $$

Here, $\rho_v$ is the vapor density, $\alpha_v$ is the vapor volume frac-

tion, and the source terms, $\dot{m}_{evap}$ and $\dot{m}_{cond}$, model evaporation and condensation rates, respectively. Additionally, $u_i$ and $x_i$ are the Cartesian velocity and directional components. When combined with a Navier-Stokes-based solver, via fluid properties, such an approach enables the modeling of cavitation.

Cavitation modeling aims to approximate the source terms in the vapor conservation equation, i.e., $\dot{m}_{evap} - \dot{m}_{cond}$, and is an area that continues to progress. The initial efforts of Merkle [1] and Kunz [2] were more focused on numerics and developed heuristic models for cavitation. Later, the Rayleigh-Plesset Equation (RPE) was explored as a physics-based method to improve cavitation models. The RPE is a model that predicts nuclei growth, collapse, and the dynamics associated with it [5, 6]. The resulting nonlinear, second-order ordinary differential equation is given as:

$$ R\ddot{R} + \frac{3}{2}\dot{R}^2 + \frac{4\nu_l\dot{R}}{R} + \frac{2S}{\rho_l R} = \frac{p_v - p}{\rho_l} + \frac{p_{G_0}}{\rho_l} \left(\frac{R_0}{R}\right)^{3\gamma} \quad (2) $$

This model describes the dynamics of an isolated bubble with a radius $R$ as it is exposed to a driving pressure, $p$. Some of the other quantities in this relation are the liquid fluid properties, i.e., kinematic viscosity ($\nu_l$), density ($\rho_l$), and saturated vapor pressure ($p_v$). In addition, the nuclei properties such as its initial radius, $R_0$, and pressure, $p_{G_0}$, as well as the surface tension, $S$. This RPE is a

<footnote n="1">¹Applied Research Laboratory, Department of Aerospace Engineering, The Pennsylvania State University</footnote>

<footnote n="2">²Department of Mechanical Engineering, The Pennsylvania State University</footnote>

<footnote>*Corresponding author: mpk176@alumni.psu.edu</footnote>
physics-based model that gives accurate predictions of
cavitation behavior.

In the context of Eulerian-based cavitation models,
reduced forms of the RPE were explored in the efforts of
Singhal et al. [4], Sauer and Schnerr [3], and Zwart et al.
[7]. These models all aim to model evaporation and con-
densation processes using a reduced RPE. One issue with
these models is that modeling evaporation and condensa-
tion processes are processes completely neglected in the
RPE. Such a modeling discrepancy still remains within
the Eulerian CFD models for cavitation. Alternatively,
Lagrangian modeling efforts of the cavitation process also
exist [8, 9] and do so without the aforementioned model-
ing mismatch. These Lagrangian efforts directly model
the cavitation processes using RPE and do in fact predict
nuclei growth and collapse. Despite their progress and
improved physical accuracy, these Lagrangian models
remain as an isolated branch of cavitation modeling from
the Eulerian efforts. Some efforts seem to be moving
to better include such cavitation dynamics as indicated
by Wang and Brennen [10] for quasi-one-dimensional
flows and Balu [11, 12] for flows without advection. In
addition, recent work [13] evaluated traditional cavita-
tion model formulations with respect to the RPE. In
the evaluations, it was found that the present model-
ing approaches are all first-order representations of the
cavitation and neglect several terms that lead to inaccu-
racy. In this work, we focus on exploring CFD cavitation
models with respect to cavitation dynamics as described
using by the RPE. The present efforts essentially push
the understanding of CFD-based cavitation modeling to
better incorporate bubble dynamics (from a growth and
transport perspectives).

A second aspect of the paper intends to focus on the
transport of the cavitation bubbles. Specifically, the ac-
curacy of the homogenous model assumption is evaluated.
This intends to address the issue of whether or not the
velocity vector in Eqn. (1) should be the same for the
bubble and liquid phases. This fundamental question,
common to boiling flows, that to the authors’ knowledge
is yet to be addressed in the context of cavitation. Com-
bining and utilizing these two approaches, we believe we
are working to develop an improved cavitation model.
The present paper is outlined as follows. The numerical
methods and test cases used for the final assessments
are presented. This is followed by a reformulation of the
cavitation model and assessments of its accuracy with
respect to the RPE for simple isolated bubble scenarios.
Next, we present a two-fluid modeling approach for a
cavitating flow with bubble slip with respect to the liquid.
This is followed by a comparison of the various model
formulations to evaluate the overall sensitivity of the
effect of cavitation modeling and bubble slip. Lastly, the
findings from these efforts are summarized with respect
to future cavitation modeling.

# 1. NUMERICAL METHODS

The present study is developed in the context of the
commercial CFD code, Star-CCM+ [14]. In this work,
we use a two-fluid, Eulerian-based multiphase-modeling
approach that aims to conserve phase mass through
a volume-of-fluid-like formulation. The CFD model
is based on a pressure-based, segregated-flow model
that conserves mass and momentum of two different
phases. For the liquid water phase, an incompressible
model is used where the liquid properties are as fol-
lows: ρ<sub>l</sub> = 1000kg/m<sup>3</sup>, μ<sub>l</sub> = 0.001Pa-s. The gas phase
assumes an incompressible gas model with properties sim-
ilar to air, i.e., ρ<sub>v</sub> = 1.2kg/m<sup>3</sup>, μ<sub>v</sub> = 1.85 × 10<sup>-5</sup>Pa-s.
Note that air properties are used as they are both similar
to vapor and, as discussed through this work, cavitation
involves mixtures of nuclei and vapor gas. Hence, the
fluid properties relevant to cavitation models is not pure
vapor and this assumption should not affect the accu-
racy of the model. Throughout this work, each phase
consists of its own velocity vector. Note that the ho-
mogenous model is achieved using drag values that drive
the slip to very low values. The numerical scheme is
formally second-order accurate in space (upwind, finite
volume) and time (backward difference). In terms of the
turbulence, a Reynolds Averaged Navier-Stokes (RANS)
turbulence model is used. The model is based on the
Spalart-Allmaras turbulence model [15]. Note that the
present analyses are independent of turbulence model
and findings should extend to any turbulence model
choice.

The numerical tests are performed using a generic,
two-dimensional, test case. For this effort, a NACA 0012
was chosen. For the case considered, the model had
a chord-based Reynolds number (Re<sub>c</sub> = ρ<sub>l</sub>cV<sub>∞</sub>/μ<sub>l</sub>), of
5 × 10<sup>6</sup>, a cavitation number (σ<sub>v</sub> = (p − p<sub>v</sub>)/(1/2ρ<sub>l</sub>V<sub>∞</sub>)),
of 2.5, and an angle of attack of 10 deg. The numerical
mesh was developed using a hex-based, two-dimensional
mesh generator internal to Star-CCM+ [14]. The mesh is
depicted in Fig. 1. In the left part of Fig. 1, the full CFD
domain is plotted. The far-field boundary conditions are
roughly 10 chords lengths from the airfoil on all four sides.
A picture of the CFD mesh near the airfoil is provided
in the right part of Fig. 1, where the prism layer can be
observed. The mesh used in this effort uses an adaptive
wall function and in this work, we use a mesh with a wall
y<sup>+</sup> (i.e., the dimensionless wall spacing) of roughly 1.0
at the trailing edge, and has roughly thirty-two prism
layers. Lastly, the mesh has roughly 300 points resolving
the airfoil. Note that in studying bubbles in the context
of a two-fluid model interacting with a wall will require
further evaluation beyond the scope of the present effort
[16]. In general, this mesh is reasonable for cavitation
and is certainly sufficient for this purpose of this study.
<figure><img src="image_3.png" /><figcaption>Figure 1. Computational mesh used in the CFD studies throughout this effort. The left figure plots the full CFD domain and the right figure plots a blowup near the hydrofoil.</figcaption></figure>

## 2. REFORMULATED CAVITATION MODEL

### 2.1 Cavitation Model Formulation

The first aspect of the present effort aims to improve cavitation modeling by more faithfully representing the RPE. In the context of a VOF-like formulation, the vapor mass conservation equation can be written as [1, 2]:

$$ \frac{\partial(\rho_{v,eff}\alpha_v)}{\partial t} + \frac{\partial(\rho_{v,eff}\alpha_v u_i)}{\partial x_i} = \rho_{v,eff}(S_{V,evap} - S_{V,cond}) \quad (3) $$

Note the source terms are now referred to as $\rho_{v,eff}(S_{V,evap} - S_{V,cond})$, where $S_{V,evap}$ and $S_{V,cond}$ respectively represent the volume growth and shrinkage of nuclei. This form is preferred as the RPE is a volume change model and does not model phase change. The effective vapor density, $\rho_{v,eff}$, is introduced as an effective density of the vapor for which the momentum is based upon. In general, the true density cannot be used as cavitation is both compressible and composed of multiple gas species. To investigate improving cavitation models, we use a process that transforms from $R$ to $\alpha_v$. Note that $R$ is the radius of an isolated bubble and is the dependent parameter in the RPE. With the use of a physically-based input, i.e., bubble concentration, $N_b$, there is a direct relation between volume fraction and bubble radius [13], which is given as

$$ \alpha_v = N_b \frac{4}{3} \pi R^3. \quad (4) $$

Note that $N_b$ is described as the number of bubbles per unit volume. Similarly, the rate-of-change in volume fraction can be described (for an incompressible model) as

$$ \dot{\alpha}_v = N_b \frac{4}{3} \pi \frac{d}{dt}(R^3) = 4\pi N_b R^2 \dot{R}. \quad (5) $$

Similarly, the second-rate-of-change in volume fraction can be described as

$$ \ddot{\alpha}_v = N_b 4\pi \frac{d}{dt}(R^2 \dot{R}) = 4\pi N_b (R^2 \ddot{R} + 2R \dot{R}^2). \quad (6) $$

Such transformations enable cavitation models to be directly based on the RPE. Although there are many choices to incorporate RPE into CFD-based cavitation models [3, 4, 7, 8], what the present effort aims to find are the key terms associated with the RPE with respect to CFD cavitation models. In the work of Kinzel et al [13], analytic solutions of bubble dynamics from Brennen [17] are used as reference solutions. In the present work, the application of those analytic solutions is extended to provide an alternative basis for the formulation of a cavitation model. Here, the radial growth rate of a bubble is estimated using an analytic solution to the RPE with the assumption of a step change in the pressure, yielding[17]:

$$ \dot{R} = \left( \frac{2}{3} \frac{(p_v - p)}{\rho_l} \left[ 1 - \frac{R_0^3}{R^3} \right] + \frac{2}{3} \frac{p_{g0}}{\rho_l} \frac{1}{1-\gamma} \left[ \frac{R_0^{3\gamma}}{R^{3\gamma}} - \frac{R_0^3}{R^3} \right] + \frac{2S}{\rho_l R} \left[ 1 - \frac{R_0^2}{R^2} \right] \right)^{1/2} \quad (7) $$

Note that in this model form, the viscosity term was neglected to arrive at this analytic solution[17]. In addition, in outside of nucleation and collapse events, $R_0 << 1$, surface tension is of lower importance, hence, the relation reduces to

$$ \dot{R} = \left( \frac{2}{3} \frac{(p_v - p)}{\rho_l} + \frac{2}{3} \frac{p_{g0}}{\rho_l} \frac{1}{1-\gamma} \left[ \frac{R_0^{3\gamma}}{R^{3\gamma}} - \frac{R_0^3}{R^3} \right] \right)^{1/2} \quad (8) $$

This analytic solution provides a RPE-based model that can help construct a physically-based cavitation model for bubble growth. In terms of the driving terms (in the root), present cavitation models are limited to the first term, i.e., $(p_v - p)/\rho_l$, and the remaining terms are neglected. When the model is converted into a source
term, the following model is obtained

$$S_{V, evap} = C_{evap} C_{\infty} \\ \times \alpha_v^{2/3} \left( \frac{2}{3} \frac{(p_v - p)}{\rho_l} + \frac{2}{3} \frac{p_{g0}}{\rho_l} \frac{1}{1 - \gamma} \left[ \frac{\alpha_0^{\gamma}}{\alpha^{\gamma}} - \frac{\alpha_{v,0}}{\alpha} \right] \right)^{1/2}, \quad (9)$$

where the nucleation-site concentration is assumed to be constant and yields its associated constant as

$$C_{\infty} = \frac{N_b^{1/3}}{6^{2/3}} \pi^{1/3}. \quad (10)$$

Note that this is a model that is distinctly different than previous works [1, 2, 3, 4, 7].

Now consider the bubble collapse events. These models are then, similarly, modeled using a model for bubble collapse [17] given as

$$\dot{R} = - \left( \frac{R_{max}}{R} \right)^{3/2} \left( \frac{2}{3} \frac{(p-p_v)}{\rho_l} + \frac{2}{3} \frac{p_{g0}}{\rho_l} \frac{1}{1-\gamma} \frac{R_0^{3\gamma}}{R^{3\gamma}} - \frac{2S}{\rho_l R} \right)^{1/2} \quad (11)$$

Neglecting surface tension, this approach reduces to

$$\dot{R} = - \left( \frac{R_{max}}{R} \right)^{3/2} \left( \frac{2}{3} \frac{(p-p_v)}{\rho_l} + \frac{2}{3} \frac{p_{g0}}{\rho_l} \frac{1}{1-\gamma} \frac{R_0^{3\gamma}}{R^{3\gamma}} \right)^{1/2} \quad (12)$$

Then converting to CFD-variables, a model for the bubble growth can be obtained as

$$S_{V,cond} = C_{cond} C_{\infty} \alpha_v^{2/3} \left( \frac{\alpha_{max}}{\alpha} \right)^{1/2} \times \left( \frac{2}{3} \frac{(p-p_v)}{\rho_l} + \frac{2}{3} \frac{p_{g0}}{\rho_l} \frac{1}{1-\gamma} \frac{\alpha_{v,0}^{\gamma}}{\alpha^{\gamma}} \right)^{1/2}, \quad (13)$$

One last issue is the $\alpha_{v,max}$ term, which is a time-history-related quantity that is not an algebraically computed value at a given point in the flow field. The challenge of this term arises as there is a time history associated with it and obtaining such a quantity is not well suited to CFD models without an auxiliary transport equation. For the current state of the model, it is assumed that $\alpha_{v,max} = 1.0$, which implies that the collapse energy is based on the maximum bubble size available. Hence, the model is expected to be accurate when processed through a sheet cavity. However, in cloud cavitation, where volume fractions do not reach unity in CFD, the model will likely over-estimate the collapse rate in cloud

cavitation. With this assumption, the model simplifies to

$$S_{V, evap} = C_{cond} C_{\infty} \alpha_v^{1/6} \times \left( \frac{2}{3} \frac{(p-p_v)}{\rho_l} + \frac{2}{3} \frac{p_{g0}}{\rho_l} \frac{1}{1-\gamma} \frac{\alpha_{v,0}^{\gamma}}{\alpha^{\gamma}} \right)^{1/2}. \quad (14)$$

The resulting model provides a CFD model that represents the RPE in the same way as the approximate models do for the bubble dynamics. Prior to proceeding, it is important to list the underlying assumptions of this modeling approach. To the best of the authors' knowledge, these assumptions are as follows:

1. The underlying assumptions of the RPE are retained (i.e., spherical bubbles, no interaction with neighboring bubbles, no thermal cavitation)

2. Cavitation process relevant to CFD behave like groups of bubbles governed by RPE

3. The model breaks down as the vapor volume fraction approaches unity (i.e., sheet cavities), and the concept of subgrid-scale isolated bubbles it no longer valid

4. Constant nuclei distribution through the fluid domain

5. Cavitation collapse goes through a process of the largest nuclei growth possible.

Keeping these assumptions in mind, one goal of this effort is to develop a more physically accurate cavitation model.

## 2.2 Cavitation Model Sensitivity

In this section, the sensitivity of various CFD cavitation models is evaluated against sensitivities in the RPE. This process is performed in the context of an ODE solution method [11, 12] to the RPE for an isolated bubble. In theory, the RPE will display sensitivities that are expected to be the most physical, hence, is used as the reference solution for evaluation of the CFD-cavitation models. Note that for these cases, $p_{g0}$ is the initial nuclei pressure and is given by $p_{\infty} - p_v + 2S/R_0$. Additionally, the viscosity is neglected, the specific gas ratio, $\gamma$, is 1.4, and the surface tension, $S$, is 0.073N/m. The test case is based on a bubble (of various radius) exposed to a pressure of 1 MPa (i.e., $p_{\infty}$). The bubble is then exposed to a sharp drop in pressure ($p_{low}$) for a duration of 300µs, then gradually returns to $p_{\infty}$ at 500µs as indicated in Fig. 2 (b). Variations (outlined in Table 1) are explored that include variable initial bubble radius and the value of the low pressure the bubble is exposed to. The results are plotted in Fig. 2. In Fig. 2 (a), are comparisons of the various models to the reference RPE equation. The constants for the Singhal model is tuned to reasonably represent the RPE for this condition. The constants from the Kunz model are based on the Singhal constants and adapted based on the rules developed in
Kinzel et al. [13]. The overall results indicate that, when the cavitation constants are correctly established, the Singhal and Kunz models replicate the RPE. Note that this is also expected to occur for the Sauer, Zwart and Merkle models. The result of Fig. 2 (a) is equivalent to establishing CFD-cavitation model constants for a single condition, and a test of the fidelity of the model is understanding how well the model extrapolates to other conditions.

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison of the predictions from the RPE, Singhal, Kunz, and Reformulated models for the growth and collapse of an isolated bubble exposed to pressure in part (b). This plot indicates that the cavitation models can replicate the RPE.</figcaption></figure>

The next aspect of evaluating the cavitation models aims to understand the sensitivity of the model to changes in the empirical constants. This is assessed for the Singhal, Kunz, and Reformulated models in Fig. 3. These evaluations are all performed for Case 1 conditions. Hence, the results from Fig. 2 (a) care co-plotted with corresponding results when the cavitation-model constants are increased and decreased by 20%. As indicated in the figure, for each cavitation model, the predictions change only a moderate amount.

<table><caption>**Table 1.** Parameters used in evaluations</caption><thead><tr><th rowspan="2">Case Name</th><th colspan="2">—Inputs—</th></tr><tr><th>R<sub>0</sub></th><th>p<sub>low</sub></th></tr></thead><tbody><tr><td>Case 1</td><td>20 μm</td><td>0.1p<sub>v</sub></td></tr><tr><td>Case 2</td><td>40 μm</td><td>0.1p<sub>v</sub></td></tr><tr><td>Case 3</td><td>20 μm</td><td>0.9p<sub>v</sub></td></tr><tr><td>Case 4</td><td>20 μm</td><td>2.0p<sub>v</sub></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Comparison of the predictions from the RPE (black-dashed line), Singhal (upper plot), Kunz (center), and Reformulated (lower) models for the growth and collapse of an isolated bubble exposed to a pressure similar to the profile in 2(a) (b). This plot highlights the sensitivity of model results to the cavitation constants for Case 1. The dotted lines indicate the results when changing the cavitation constants by ±20%.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Comparison of the predictions from the RPE, Singhal, Kunz, and Reformulated models for the growth and collapse of an isolated bubble exposed to a pressure similar to the profile in 2(a) (b). This plot indicates the impact of nuclei size (Case 2, R<sub>0</sub> = 40μm, p<sub>low</sub> = 0.1p<sub>v</sub>).</figcaption></figure>

In Figs. 4-6, are evaluations of the cavitation constants for various conditions. In Fig. 4, the same pressure (as in Case 1) is applied, however, the nuclei size is increased from an initial radius of 20μm to 40μm. Such a condition may represent a sensitivity to water quality such as water tunnel conditions as compared to seawater. The result of the plot clearly indicates that the Singhal and Kunz models have no sensitivity to this physical characteristic. However, the RPE and reformulated models both display sensitivities. These studies are further evaluated with respect to the pressure. In Fig. 5 is a condition where the models are applied to condi-
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Comparison of the predictions from the RPE, Singhal, and Kunz models for the growth and collapse of an isolated bubble exposed to a pressure similar to the profile in 2(a) (b). This plot indicates the impact of the minimum pressure (Case 3, $R_0 = 20\mu m$, $p_{low} = 0.9p_v$).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison of the predictions from the RPE, Singhal, and Kunz models for the growth and collapse of an isolated bubble exposed to a pressure similar to the profile in 2(a) (b). This plot indicates the impact of nuclei exposed to a minimum pressure above vapor pressure (Case 4, $R_0 = 20\mu m$, $p_{low} = 2.0p_v$).</figcaption></figure>

tions where the low-pressure event is only slightly below vapor pressure. As a cavitation model tends to drive cavitation based on a fraction of the cavitation number, i.e., $\sigma_v = (p_v - p)/(1/2\rho_l V_\infty^2)$, hence, the higher the speed, the lower the pressure gets. Such a sensitivity may correlate to a lower speed cavitation. The RPE model displays only a slight decrease in the resulting radius (when compared to Fig. 2 (a)), whereas the Singhal and Kunz models tend to predict a much smaller cavitation event, hence, fail to predict this sensitivity. The reformulated model essentially reproduces the sensitivity of the RPE for this condition. In Fig. 6, is a final comparison that evaluates nuclei cavitating at pressures higher than vapor pressure. In this case, the RPE clearly shows a cavitation will occur. However, the Singhal, Kunz, and reformulated models fail this test. Hence, despite improving the state of CFD-cavitation modeling, cavitation occurring above vapor pressure is not captured. Nevertheless, from an individual bubble perspective, the reformulated cavitation improves the state of the art in terms of replicating the RPE behavior.

An overview of the sensitivities evaluated is provided in Fig. 7. This plot indicates the error in the resulting maximum diameter for each model, with respect to that predicted by the RPE. Recall that error in the first case

indicates the error in the established model constants. The error in the remaining cases indicates an error in the ability to predict the same sensitivities as does RPE. In this plot, and for this test problem, it is clear the Singhal and Kunz models, both, require modification of the cavitation constants to capture these effects. The reformulated model appears to capture the sensitivities quite well. And, none of the models capture the fourth case, where the low pressure remains above vapor pressure.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Quantified comparisons of the predictions from the RPE, Singhal, and Kunz models for maximum bubble radius based on the results from Fig. 4-6. Note that the characteristics of the bubble are altered for each plot and the Singhal- and Kunz-model constants are tuned for Case 1.</figcaption></figure>

## 2.3 Modeling Cavitation-Bubble/Liquid Slip

In the context of developed cavitation, most CFD is carried out using Eulerian continuum mixture (VOF-like) approaches. The second goal of this effort aims to directly evaluate homogeneous Eulerian approaches, where local mixture equilibrium assumptions are assumed (i.e., $u_{i,v} = u_{i,l}$). One evaluation approach to assumptions involves testing model sensitivities to including more physically accurate terms. In this work, a Eulerian multiphase model formulation of the unsteady, incompressible, RANS equations as developed and implemented in Star-CCM+[14] are explored. The resulting multiphase model yields a unique velocity vector, $\bar{V}_q$, for each phase, q, at each computational point. Such a model bases mass conservation (in Cartesian-tensor notation) for each phase as

$$ \frac{\partial(\rho_q \alpha_q)}{\partial t} + \frac{\partial(\rho_q \alpha_q u_{i,q})}{\partial x_i} = S_q \quad (15) $$

with $\rho_q$, $u_{i,q}$, and $\alpha_q$ representing the density, velocity, and local volume fraction, respectively, of each phase q. Here, $S_q$ is the source that models cavitation, (i.e., $S_q = \rho_v (S_{V,evap} - S_{V,cond})$). Note that for N-phases, the following volume conservation constraint is satisfied:

$$ \sum_{q=1}^{N} \alpha_q = 1.0 \quad (16) $$
Momentum conservation for each phase (both gas and liquid) is expressed as

$$ \frac{\partial(\rho_q \alpha_q u_{j,q})}{\partial t} + \frac{\partial(\rho_q \alpha_q u_{i,q} u_{j,q})}{\partial x_i} = -\alpha_q \frac{\partial p}{\partial x_j} + \frac{\partial \tau_{i,j,q}}{\partial x_j} + F_{Hyd,j,q,p} \quad (17) $$

In this relation, $F_{Hyd,j,q,p}$ is the hydrodynamic-based momentum exchange occurring on phase q through hydrodynamic forces with the other phase, p. A realistic two-fluid model that captures slip demands a reasonable model of the drag behavior.

Formulation of the drag force typically requires a drag model that is based on a normalized drag coefficient. The hydrodynamic force for such a model can be described as

$$ F_{Hyd,j,q,p} = K_{j,q,p}[u_{j,q} - u_{j,p}] \quad (18) $$

The model for the momentum exchange coefficient, $K_{j,q,p}$, is defined according to the content of vapor bubbles. The model form is described as

$$ K_{j,q,p} = 0.75 \frac{\alpha_q \alpha_p \mu_l}{d_v} Re_{q,p} c_D, \quad (19) $$

where the drag coefficient is given by the modified Schiller-Naumann [18] model given as

$$ c_D = k_B \frac{24}{Re_{q,p}} (1.0 + 0.15 Re_{q,p}^{0.687}). \quad (20) $$

Here, $Re_{q,p}$ is the Reynolds number based on the carrier-fluid properties, slip velocity, and the size of the bubble, and is defined as

$$ Re_{q,p} = \frac{\rho_l d_v V_{slip}}{\mu_l}, \quad (21) $$

where $V_{slip} = |V_v - V_L|$. Also, $k_B$ is a modification constant that is applied to account for bubbles. Lastly, the bubble diameter is the last unclosed parameter, which is estimated using spherical assumptions and is given as

$$ d_v = 2 \sqrt[3]{\frac{3 \alpha_v}{4\pi N_B}} \quad (22) $$

Such a model is fully closed and provides a platform to evaluate the impact of the homogenous model for cavitation prediction.

Before moving forward, it is worth a discussion of the two-fluid model drag model with respect to cavitation. In general, the larger value of $c_D$, and lighter the fluid enables drives the model to conditions that are homogeneous. In this work, and for consistency, a homogenous model is obtained by setting the bubble diameter to $1\mu m$, essentially providing a large $c_D$ value that forces velocity equilibrium. With respect to the model that

aims to capture slip, a model for the bubble drag is constructed. Presently available bubble drag models are not clearly applicable to cavitation [19, 20]. Specifically, the main dependency, gravity, is not directly applicable to cavitation. Hence, the current model modified the Schiller-Naumann model using a value of $k_B$ of 16/24, which is based on the efforts of Kay and Nedderman [21] of the drag behavior in the Stokes regime. This model is reasonable for this initial analysis but certainly, needs advancement moving forward.

Using the formulated two-fluid model, it is quite easy to expose the inherent weaknesses of a homogenous model for developed cavitation. For example, the simulation results plotted in Fig. 8 highlight the importance of slip using a two-fluid model (b) and (c) compared to a homogeneous mixture model (a). In these results for a 2D developed cavitation over the NACA 0012 hydrofoil case previously discussed are examined with the Singhal [4] mass transfer model. Even for this simple model (single bubble field, drag only), significant differences are observed between the two, with normalized slip velocities reaching 70% adjacent to the re-entrant jet, and significantly more unsteady richness and cavity diffusion arising in the conventional homogenous model. Specifically, the driver of this slip is found to relate to high curvature and decreased interactional forces between the bubbles and liquid. As indicated in Fig. 9, as the nuclei show to grow in size in part (a), they have a corresponding increase in Reynolds number (b) and decrease in drag coefficient (c). This scenario leads to decreased forcing on the bubbles, while simultaneously increasing the content. In Fig. 9 (d), are trajectories of the bubbles (black streamlines) as compared to the liquid (white streamlines). In comparing the streamlines, what we observe is that the bubbles reach the cavity terminus at a point when they have a corresponding low drag. Hence, the high curvature region leads to a scenario where the bubbles tend to have a high slip with respect to the liquid. This slip allows the bubbles to recirculate back into the cavity, leading to a more gas-rich cavity. In the homogenous model, where the bubbles are not allowed to have such slip, the liquid more effectively entrains the gas from the cavity. Hence, despite the simplicity of this model, and the obvious implications of neglecting vapor bubble dynamics, the homogenous assumption appears to an assumption that should be reconsidered.

# 3. DISCUSSION

## 3.1 Overall Impact of Model Reformulations

An overall comparison of the evaluated physical modeling changes is displayed in Fig. 10. In Fig. 10 is an evaluation of the impact of the cavitation model (Singhal with respect to the reformulated) and the slippage of the cavitation bubbles. Note that these are all time-averaged plots of the vapor volume fraction. Fig. 10 (a) indicates the result from the Singhal cavitation model with a
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of (a) homogeneous and (b) two-fluid cavitation predictions for a NACA 0012, $\alpha = 10$ deg, $Re_c = 5 \times 10^6$, $\sigma_v = 2.4$. For (a) and (b), the contour lines plotted are pressure coefficient ($C_p$) and filled contours are the gas volume fraction ($\alpha_v$) at three corresponding timesteps. In part (c) are normalized slip velocity for the two-fluid model, at a given timestep.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Evaluation of the bubble-liquid slippage. Results indicate an increase in bubble size and Re, hence, drag coefficient drops. The trajectories of the bubbles (black streamlines) as compared to the liquid (white streamlines) indicates slippage is strongest in high curvature regions near the aft end of the cavity.</figcaption></figure>
homogenous model. When compared to Fig. 10 (b), which has the Singhal model with the two-fluid model, the impact of bubble slippage is observed. Note that the slippage leads to a tendency for the vapor to be sustained in the cavity leading to an increased vapor volume fraction. In comparing Fig. 10 (a) and (c), where (c) has the reformulated-homogenous model, presents the impact of the reformulated cavitation model. This model is indicating a moderate impact to the cavity character, however, as observed in Figs. 2 to 5, the impact of the reformulated cavitation model is most felt in the extending to other flow conditions (which is not directly evaluated in this plot). Finally, comparison of Fig. 10 (c) and (d) indicates the impact of the two-fluid modeling with the reformulated cavitation model. In this case, a dramatic impact to the predicted cavity is observed. The cause of this is not yet fully understood, but it likely relates to the reformulated cavitation model predicting larger bubbles with more slip that lead to additional dispersion of the cavity. The overall result of this study clearly suggests that cavitation modeling has several first-order sensitivities to improved physical modeling that require additional evaluations.

## 3.2 Cavitation Model Assumptions

At this point, it is worth revisiting the basic RPE assumptions with respect to the results for a cavitating flow around a hydrofoil. When cavitation occurs at high void fractions (i.e., $\alpha_v > 0.05$), many of the underlying assumptions break down. From Fig. 10, this scenario is clearly the case throughout the cavity. In this such high-gas-content regions, we fully expect (1) breakdown of isolated bubbles, (2) non-spherically symmetric bubbles, and (3) breakdown of the constant nuclei concentration. For these reasons, the accuracy of the reformulated model requires caution. From a physical viewpoint, the model form for groups of bubbles described in Brennen [17] has a similar model form as the RPE. In addition, it is expected that a non-spherical bubble will respond much like a spherical bubble. Hence, cavitation constants may be a suitable approach to correct such model discrepancies. Of more concern is the breakdown of the constant nuclei concentration. Specifically, when the cavity grows, the nuclei merge. In addition, during the collapse of a cavity, the cavity break-down process is very likely to modify the nuclei concentration (i.e., nuclei content is expected to increase through cavity-collapse events). Such non-linear events are not expected to be corrected with an empirical constant. Lastly, at some point in a cavitating flow, vapor extraction is expected to be of importance. Especially in stationary cavitation events where the cavity has time to generate vapor. Hence, there are a number of fundamental modeling issues in cavitation model that remain.

## 3.3 Perspectives on Cavitation Modeling

In these efforts, several modified perspectives relevant to CFD cavitation modeling are developed. The first perspective pertains to the common conception of mass transfer processes occurring in cavitation. As suggested in Brennen [17], cavitation is dominated by gas expansion rather than mass transfer. For this reason, the cavitation modeling efforts should physically correspond to this and should consider departing from being modeled as boiling. Hence, rather than referring to the processes as mass transfer, we prefer to discuss them as volumetric source terms as follows:

$$ \frac{\partial(\rho_g^* \alpha_g^*)}{\partial t} + \frac{\partial(\rho_g^* \alpha_g^* u_i)}{\partial x_i} = \rho_g^*(S_{V, evap} - S_{V, cond}) \quad (23) $$

Here, $\rho_g^*$ is an effective density of the gas/vapor mixture that the model sees and need not correspond to the actual gas/vapor density (but should remain similar for momentum purposes). The modified volume fraction, $\alpha_g^*$, should also represent the gas/vapor mixture and not just a mixture. Such a model is more consistent with both cavitation and basic RPE cavitation model formulation (herein as well as other cavitation models [1, 2, 3, 4, 7]).

The aspect of modeling cavitation that needs to be considered is model formulation with respect to compressible and incompressible flows. Previous works (e.g. Kinzel et al. [22] and others) tend to use cavitation models developed and validated for incompressible flow. The reformulated model presented herein clearly indicates the importance of compressible terms. In the context of an incompressible code, these compressible terms are reasonably modeled. However, in the context of a compressible fluid flow, these terms certainly need to be removed (or partially removed). Such evaluations are critical to the future success of modeling compressible cavitation.

Lastly, cavitation with thermodynamic effects is also another concern. As stated throughout this effort, cavitation bubbles grow more from gas expansion than they do from phase change. When modeling thermodynamic effects, such as done by Kinzel et al. [22] and others, such factors are not accounted for. Hence, the modeling of thermodynamic effect with respect to cavitation needs to be revisited.

## 4. CONCLUSIONS

The present effort revisited cavitation modeling with respect to the RPE. Rather than reducing the RPE, the present approach directly utilized analytic solutions to the RPE. The efforts found several relevant terms tend to be missing in most CFD cavitation models. The effort displayed that the missed physics (in traditional CFD cavitation models) did not display the same physical sensitivities than those inherent to the RPE displayed. The reformulated cavitation model that considered these effects displayed the ability to better capture the RPE
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison of the homogeneous model formulations (left column) to the two-fluid models that account for bubble slippage (right column). This process is evaluated for both the Singhal model (top row) and the reformulated cavitation model (bottom row).</figcaption></figure>

sensitivities. Hence, the model may provide a more increased fidelity to CFD cavitation models.

In addition to exploring cavitation model development, the concept of including the effect of slip between the cavitation bubbles and the liquid was evaluated. The present results clearly indicated a strong impact of slip to the transport, hence, character of the cavity. Specifically, it was observed that the cavitation bubbles tended to better sustain within the cavity. Such a model sensitivity indicates that bubble/liquid slip is likely a factor demanding large computational meshes for cavitation.

In addition to these evaluations, the overall effort of the paper yield several insights. The first insight is that cavitation modeling in CFD needs to reconsider the representation of the cavitation models. Cavitation is more of a process of gas expansion/contraction than evaporation/condensation, hence, referring as these terms as mass transfer terms are not consistent with the physical phenomenon. Second, these source terms are dominated by volumetric changes and are more of an underlying representation of underlying compressible effects. Because of this, special attention needs to be taken for compressible CFD modeling of cavitation. Lastly, the present state of thermal cavitation, especially, is exposed to issues when referring to volumetric growth as a mass transfer.

## ACKNOWLEDGMENTS

This is unfunded work performed at The Pennsylvania State University, Applied Research Laboratory.

## REFERENCES

[1] Charles L Merkle, J Feng, and Phillip E O Buelow. Computational modeling of the dynamics of sheet cavitation. In *3rd International symposium on cavitation*, Grenoble, France, volume 2, pages 47–54, 1998.

[2] Robert F. Kunz, David A. Boger, David R. Stinebring, Thomas S. Chyczewski, Jules W. Lindau, Howard J. Gibeling, Sankaran Venkateswaran, and T. R. Govindan. Preconditioned Navier-Stokes method for two-phase flows with application to cavitation prediction. *Computers & Fluids*, 29(8):849–875, 2000.

[3] Jürgen Sauer and G. H. Schnerr. Development of a new cavitation model based on bubble dynamics. *Zeitschrift für Angewandte Mathematik und Mechanik*, 81(S3):561–562, 2001.

[4] Ashok K Singhal, Mahesh M Athavale, Huiying Li, and Yu Jiang. Mathematical Basis and Validation of the Full Cavitation Model. *Journal of Fluids Engineering*, 124(3):617–624, aug 2002.

[5] Lord Rayleigh. VIII. On the pressure developed in a liquid during the collapse of a spherical cavity. *Philosophical Magazine Series 6*, 34(200):94–98, 1917.

[6] Milton S Plesset. The dynamics of cavitation bubbles. *Journal of Applied Mechanics*, 16(3):277–282, 1949.

[7] Philip J Zwart, Andrew G Gerber, and Thabet Belamri. A Two-Phase Flow Model for Predicting Cav-
itation Dynamics. In *ICMF 2004 International Conference on Multiphase Flow*, 2004.

[8] Jingsen Ma, Chao-Tsung Hsiao, and Georges L Chahine. Modelling Cavitating Flows using an Eulerian-Lagrangian Approach and a Nucleation Model. *Journal of Physics: Conference Series*, 656:012160, 2015.

[9] Chao-Tsung Hsiao, Jingsen Ma, and Georges L Chahine. Multi-Scale Two-Phase Flow Modeling of Sheet and Cloud Cavitation. *30th Symposium on Naval Hydrodynamics*, 2014.

[10] Yi-Chun Wang and C E Brennen. One-Dimensional Bubbly Cavitating Flows Through a Converging-Diverging Nozzle. *Journal of Fluids Engineering*, 120(1):166-170, mar 1998.

[11] Asish Balu. *A Singularity Handling Approach to The Rayleigh Plesset Equation*. PhD thesis, The Pennsylvania State University.

[12] Asish Balu, Michael P Kinzel, and Scott T Miller. A Singularity Handling Approach for the Rayleigh-Plesset Equation. *International Journal of Numerical Methods and Applications*, 15(2):151, 2016.

[13] Michael P Kinzel, Jules W Lindau, and Robert Kunz. A Unified Model for Cavitation. *FEDSM17 2017 ASME*.

[14] CD-Adapco. Star-CCM+ 9.06, 2014.

[15] Philippe R. Spalart and S R Allmaras. A one-equation turbulence model for aerodynamic flows. In *30th aerospace sciences meeting and exhibit*, page 439, 1992.

[16] Robert F Kunz, Howard J Gibeling, Martin R Maxey, Gretar Tryggvason, Arnold A Fontaine, Howard L Petrie, and Steven L Ceccio. Validation of Two-Fluid Eulerian CFD Modeling for Microbubble Drag Reduction Across a Wide Range of Reynolds Numbers. *Journal of Fluids Engineering*, 129(1):66-79, may 2006.

[17] E. Christopher, Brennen. *Cavitation and bubble dynamics*. Oxford University Press, 1995.

[18] L Schiller and A Naumann. A drag coefficient correlation. Vdi Zeitung, 77: 318-320, 1935. *Cit{é} page*, 38.

[19] I. Roghair, Y.M. Lau, N.G. Deen, H.M. Slagter, M.W. Baltussen, M. Van Sint Annaland, and J.A.M. Kuipers. On the drag force of bubbles in bubble swarms at intermediate and high Reynolds numbers. *Chemical Engineering Science*, 66(14):3204-3211, jul 2011.

[20] M. Simonnet, C. Gentric, E. Olmos, and N. Miedoux. CFD simulation of the flow field in a bubble column reactor: Importance of the drag force formulation to describe regime transitions. *Chemical

*Engineering and Processing: Process Intensification*, 47(9-10):1726-1737, sep 2008.

[21] John Menzies Kay and Ronald Midgley Nedderman. *Fluid mechanics and transfer processes*. CUP Archive, 1985.

[22] Michael P Kinzel, Jules W Lindau, and Robert F Kunz. An examination of thermal modeling affects on the numerical prediction of large-scale cavitating fluid flows. In *7th International symposium on cavitation (CAV2009)*, Ann Arbor, USA, 2009.