# A robust and accurate finite element framework for cavitating flows with
fluid-structure interaction

Suraj R. Kashyap<EMAIL_ADDRESS>Rajeev K. Jaiman<EMAIL_ADDRESS>Department of Mechanical Engineering, The University of British Columbia,
Vancouver, BC V6T 1Z4

###### Abstract

In the current work, we present a unified variational mechanics framework for
the cavitating turbulent flow and the structural motion via a stabilized
finite element formulation. To model the finite mass transfer rate in
cavitation phenomena, we employ the homogenous mixture-based approach via
phenomenological scalar transport differential equations given by the linear
and nonlinear mass transfer functions. Stable linearizations of the finite
mass transfer terms for the mass continuity equation and the reaction term of
the scalar transport equations are derived for the robust and accurate
implementation. The linearized matrices for the cavitation equation are
imparted a positivity-preserving property to address numerical oscillations
arising from high-density gradients typical of two-phase cavitating flows. The
proposed formulation is strongly coupled in a partitioned manner with an
incompressible 3D Navier-Stokes finite element solver, and the unsteady
problem is advanced in time using a fully-implicit generalized-$\alpha$ time
integration scheme. We first verify the implementation on the benchmark case
of Rayleigh bubble collapse. We demonstrate the accuracy and convergence of
the cavitation solver by comparing the numerical solutions with the analytical
solutions of the Rayleigh-Plesset equation for bubble dynamics. We find our
solver to be robust for large time steps and the absence of spurious
oscillations/spikes in the pressure field. The cavitating flow solver is
coupled with a hybrid URANS-LES turbulence model with a turbulence viscosity
corrected for the presence of vapor. We validate the coupled solver for a very
high Reynolds number turbulent cavitating flow over a NACA0012 hydrofoil
section. Finally, the proposed method is solved in an Arbitrary Lagrangian-
Eulerian framework to study turbulent cavitating flow over a pitching
hydrofoil section and the coupled FSI results are explored for the
characteristic features of cavitating flows such as re-entrant jet and
periodic cavity shedding.

###### keywords:

Cavitation , Fluid-Structure Interaction , Homogeneous-Mixture , Stabilized
Finite Element , Partitioned iterative , Pitching Hydrofoil

††journal: Computers & Mathematics with Applications

##  1 Introduction

Cavitation is ubiquitous in natural and industrial industrial systems such as
hydrofoils, nozzles, pumps, underwater vehicles and marine propellers. While
cavitation can be a major source of noise, vibration and material erosion in
these systems as unwanted effects [10, 34], useful applications of cavitation
have been developed in underwater cleaning [55, 12], ultrasonic cleaning [11,
60], biomedical procedures such as lithotripsy [5, 43] and enhanced drug
delivery [56, 27], etc. The phenomenon of cavitation involves the phase change
of liquid into vapor and a highly complex interaction between the vapor and
the liquid phases. Most liquids have inherent points of weaknesses in the form
of entrained microscopic gas bubbles, suspended particles, crevices along the
shared boundaries with solid structures, and ephemeral voids created by the
thermal motions of the liquid. When the flowing liquid encounters a region of
low pressure, it has a propensity to rupture at these locations of weaknesses
due to the tensile forces [7]. This rupturing creates cavities in the fluid,
which can be filled with the original entrained gas or by vapor generated by
evaporation at the cavity interface. These cavities can then be convected by
the flow until they encounter a region of high pressure, where they can
undergo violent collapse. Cavitation is often encountered in marine propeller
operations where fluid acceleration over propeller blades generates low-
pressure regions near the blade surface. At these locations, nuclei present in
seawater are prone to rupturing and forming vapor/gas-filled cavities. The
cavities can then remain attached to or be shed from the blade surface in the
form of cavitating structures with varying energy content [45]. The onset of
cavitation influences the fluid-structure dynamics of propeller operation and
can have several undesirable effects including performance degradation,
vibration, material erosion and noise emission [7]. Traditional approaches to
model the coupled fluid-structure dynamics in cavitating flows have often
focused on a one-way coupling between a flow solver for the cavitation
hydrodynamics and a separate finite element solver for structural deformation.

The reduction of noise emission in marine vessels is of interest both from an
industrial and a marine-environmental perspective. For example, [58] showed
that in the cavitating regime, propeller noise dominates all other sources of
self-noise from ships, including electrical noise, machinery noise and
boundary layer noise. Recently, [10] provided an excellent review of noise
from cavitating propellers and identified two broad categories: (i) a
broadband noise component resulting from the sudden collapse of cavities and
vortices, and (ii) tonal noise components from periodic fluctuations in the
cavity volumes. In a classical work, [34] identified that tonal noise emission
in propellers was a result of blade vibration due to irregular cavitation and
vortex-shedding dynamics. Several experimental studies [10] have since
confirmed this. In addition, [3, 39, 4] showed that resonance in cavity-filled
vortices shed from the blade tip can also emit intense tonal noise. Thus, a
numerical study of propeller cavitation noise needs to consider this complex
hydrodynamic interplay between cavitation and vortex shedding, and their FSI
effects with the propeller blades. Figure 1 demonstrates some of these
predominant noise-generating mechanisms in hydrodynamic cavitation, with
$\Omega^{\mathrm{f}}$ representing the fluid domain, $\Omega^{\mathrm{s}}$ the
solid structural domain and $\Gamma^{\mathrm{fs}}$ the fluid-structure
interface.

$\Omega^{\mathrm{f}}$$\Omega^{\mathrm{s}}$$\Gamma^{\mathrm{fs}}$collapsing
cavitiesresonating tip vortex cavityre-entrant jetshed cavitiesleading-edge
partial cavities$\&$ vortices$U_{\infty}$ Figure 1: Representative schematic
demonstrating some of the prominent noise sources from cavitating blades.
$\Omega^{\mathrm{f}}$ and $\Omega^{\mathrm{s}}$ are the fluid and solid
domains, $\Gamma^{\mathrm{fs}}$ represents the interface between them. The
emitted noise can be decomposed into two main components: (i) tonal noise
originating from blade vibration and resonance in tip vortex cavities, and
(ii) broadband noise resulting from violently collapsing cavities and
vortices.

###  1.1 Transport-equation based modeling of cavitating flows

Cavitation in propellers manifests in the form of cavitating structures that
exist across multiple orders of spatial and temporal scales [4, 7], making the
study of marine cavitation a challenging task. Thus, each of the different
approaches developed for the numerical modeling of cavitating flows is
generally computationally feasible only for a select set of flow
configurations. A popular approach is to represent the fluid as a homogeneous
mixture of liquid and vapor. The homogeneous mixture-based approaches differ
in the estimation of the density field and can be broadly classified into two
categories.

The first category assumes equilibrium flow theory and the density is
calculated using equations of state (EoS) [48, 52]. For isothermal flows, a
barotropic equation of state is used to represent the density field as a
function of the pressure [16, 13]. The equilibrium flow approach has the
advantage of easier implementation. It also does not require the use of
empirical coefficients for modeling and relies on well-established equations
of state. However, these models are generally used with the compressible Euler
equations and solved using density-based solvers. In [21], the authors
compared a compressible density-based method and an incompressible pressure-
correction method to study cavitating flows with the barotropic EoS and
reported better correlations with the compressible approach. Thus, this
approach requires very small time-steps to capture the pressure wave
propagations in the compressible fluid [19]. Another limitation of the
equilibrium flow models based on barotropic EoS is that the gradients of the
density and the pressure fields are parallel. Thus, the baroclinic torque
which is proportional to the cross product $\left(\nabla\rho\times\nabla
p\right)$ of these gradients is zero. This can result in inaccurate estimates
of the vorticity production, which is an important feature of cavitating flows
particularly in the closure region of attached cavities [22]. Hence this
approach may not be appropriate for modeling the physics of cavitating flows
over hydrofoils.

The second category of homogeneous mixture-based approaches assumes the pure
liquid and vapor phases to be incompressible. The two-phase mixture density is
interpolated based on a phase indicator. This phase indicator is generally in
the form of the local phase fraction of either the liquid or the vapor phase.
The phase indicator in the computational domain is obtained as the solution of
a scalar transport equation. We refer to these transport-equation based models
as TEMs in the rest of the article. The TEMs employ a source term, which is
indicative of a finite mass transfer rate between the two phases by the
process of cavitation. TEMs generally vary in the formulation of the source
term based on phenomenological arguments. Merkle et al. [41] used dimensional
arguments for bubble clusters to relate the source term to the local pressure
and phase fraction of liquid. Kunz et al. [35] considered a similar source
term as [41] to model the evaporation process, but modified the condensation
term employing a simplified form of the Ginzburg-Landau potential. Later
models by Schnerr-Sauer [47], Zwart et al. [61] and Singhal et al. [54]
assumed the cavities to be present in the form of clusters of spherical
bubbles and used directly a simplified form of the Rayleigh-Plesset equation
[7] for spherical bubble dynamics to model the mass transfer rate. These
models differ in multiplier terms that were derived using different
phenomenological arguments for the underling bubble-bubble interaction. TEMs
have been applied to the study of several cavitating flow configurations,
including hydrodynamic cavitation over hydrofoils [49, 20, 29].

Despite wide applicability, one limitation of the TEM approach is the use of
semi-empirical coefficients in the source terms that need to be tuned for
particular flow configurations. In addition, since the density is interpolated
using the phase indicator, large spatial gradients in the density field exist
across the cavity interface. This can result in unphysical oscillations/spikes
in the pressure in the vicinity of the interface. For incompressible flow
simulations, these pressure spikes can propagate rapidly throughout the
computational domain, leading to numerical instability. These numerical
artifacts have been reported in several works such as [19], and special
treatments are required for discontinuity capturing across the interface. [50]
suggested the ability to handle these spurious pressure spikes as one of the
conditions for a robust solver for cavitating flows.

###  1.2 Review of numerical studies of two-phase FSI in marine propellers

The last decade has seen advances in the numerical study of FSI effects in
non-cavitating marine propeller operations. Simpler approaches have used
methods based on inviscid potential flow theories. [37] and [40] used a
coupled potential theory-based boundary element method (BEM) and finite
element method (FEM) approach for studying the hydro-elastic response of
flexible marine propellers. [30] presented an optimization methodology for
propellers considering FSI of the fluid and propeller blades using the panel
method (PM), which is derived from BEM. A loosely coupled PM-FEM approach is
used for the fluid-structure interaction between the hull wake and the
propeller. While inviscid models have been demonstrated to be effective for
making general design decisions for propellers because of low computational
cost, they are unable to capture the vortex-shedding process which is an
important feature of cavitating flows. In addition, multiple modeling
decisions have to be made on a case-by-case basis for different propeller
geometry and inflow conditions. Within the purview of viscous flow modeling,
[36] used a tightly coupled CFD-FEM solver to study high Reynolds number flow
over a flexible blade undergoing vortex-induced vibration. The authors
presented the requirement of tightly-coupled FSI solvers for studying large-
amplitude 3D vibrations at high Reynolds numbers.

Although hydrodynamic studies of cavitating flows over propellers abound in
literature, only limited examples of studies that consider FSI have been
found. [25] used a one-way coupled commercial computational fluid dynamics
(CFD) solver to determine hydrodynamic loads on a cavitating hydrofoil
undergoing prescribed pitching motion. [2] used a loose hybrid coupling to
couple a commercial 2D URANS solver with a 2DOF hydrofoil model to study the
effect of cavitation on the hydroelastic stability of hydrofoils. [59] used a
similar approach to study the cavity shedding dynamics and flow-induced
vibration over a hydrofoil section. Reasonable agreements were reported with
experimental measurements. However, there is a need for a robust and accurate
unified framework for the strongly-coupled FSI studies of cavitating flows
over propellers.

The finite element method lends suitably to this purpose of modeling the
caviating flows with fluid-structure interaction. It has long been staple for
the study of structural deformations, and has also been successfully applied
to fluid flow studies [32, 33]. However, not much work has yet been done on
the modeling of cavitation using the finite element method. [6] presented a
variational framework for cavitating flows applied to prescribed motions of
hydrokinetic turbines. A monolithic approach was taken for the coupling of the
cavitation TEM and flow equations. However, not much discussion has been made
on the cavity collapse pressures, shedding dynamics or the numerical stability
at high liquid-vapor density ratios seen in marine flows. Recently, Joshi and
Jaiman [32] presented a variational finite element framework to study the 3D
hydroelastic response of marine riser in turbulent marine flows, undergoing
high-amplitude vortex-induced vibrations. A strongly-coupled partitioned
approach was taken to solve the governing equations for the fluid flow, the
structural deformation and the transport of the eddy viscosity. Further, Joshi
and Jaiman [33] presented a positivity preserving variational (PPV) scheme for
the numerical solution of two-phase flow of immiscible fluids. A discrete
upwind operator was used locally in the interface region demonstrating
oscillations. This acted in the form of added diffusion, ensuring positivity
of the underlying element-level matrices. The scheme was demonstrated to be
effective in reducing spurious pressure oscillations in fluids with high-
density ratios. Numerical solutions of cavitating flows using TEMs can benefit
from similar treatment. The current work extends this framework by introducing
the physics involved in the modeling of cavitating flows.

###  1.3 Current work and contributions

In the current work, we propose a novel finite element formulation for the
numerical modeling of cavitating flows. The objective is to integrate
numerical modeling of cavitation into the framework of the stabilized finite
element methods. Two cavitation TEMs based on homogeneous flow theory are used
to model finite mass transfer rates between the liquid and vapor phases by the
process of cavitation. The cavitation TEMs used are by Merkle et al. [41] and
Schnerr-Sauer [47]. We present stable linearizations of the two models for our
proposed variational formulation. The linearized formulations are implemented
in a variational finite element framework and the elemental matrices are
imparted a positivity-preserving property for numerical stability in regions
dominated by large density ratios. The fluid flow and cavitation solvers are
coupled in a staggered partitioned manner for versatility, and predictor-
corrector iterations are used for convergence stability and accuracy. A fully-
implicit generalized-$\alpha$ time integration scheme [28] with user-
controlled high-frequency damping is used to advance the solution in time,
allowing for numerical stability at relatively coarse spatial and temporal
discretizations. Consistent with the suggestions in [50], the efficacy of the
proposed method is demonstrated with respect to the requirements for a robust
computational method for cavitating flows. The requirements can be summarized
as (i) the accurate prediction of the pressure field, (ii) the ability to
handle large density ratios of the order of 100-1000, and (iii) absence of
spurious pressure spikes across the cavity interface.

In the sections that follow, the governing equations and the proposed
formulation are first presented. Next, a benchmark test of spherical vaporous
bubble collapse is used to verify the numerical implementation. The
implementation is then validated on a turbulent cavitating flow configuration
over a hydrofoil section. Before concluding, the case of a pitching hydrofoil
is considered to explore the compatibility of the formulation with FSI and its
ability to predict flow features characteristic of cavitating flows. Finally,
we close our paper with some concluding remarks.

##  2 Numerical formulation and implementation

In this section, we present a positivity preserving variational finite element
implementation for the numerical solution of TEMs applied to cavitation
modeling. In particular, we consider two types of TEM-based homogeneous
mixture models classified as follows:

  * •

Model A: Cavitation model with nonlinear mass transfer rate by Schnerr and
Sauer [47]

  * •

Model B: Cavitation model with linear mass transfer rate by Merkle et al. [41]

In the rest of manuscript, we shall refer to the models as A and B
respectively. We design stable linearizations of models A and B for the
numerical implementation within the proposed variational finite element
framework.

###  2.1 Governing equations

The strong forms of the governing equations are presented before introducing
the variational formulation. We consider the fluid physical domain
$\Omega^{\mathrm{f}}(\boldsymbol{x}^{\mathrm{f}},t)$ with an associated fluid
boundary $\Gamma^{\mathrm{f}}(t)$, where $\boldsymbol{x}^{\mathrm{f}}$ and $t$
represent the spatial and temporal coordinates. The working fluid, consisting
of the liquid and vapor phases, is assumed to be present in the form of a
continuous homogeneous mixture. The phase indicator
$\phi^{\mathrm{f}}(\boldsymbol{x}^{\mathrm{f}},t)$ is used to represent the
phase fraction of the liquid phase at any coordinate
$(\boldsymbol{x}^{\mathrm{f}},t)$ in the homogeneous two-phase liquid-vapor
mixture. The fluid density ($\rho^{\mathrm{f}}$) and dynamic viscosity
($\mu^{\mathrm{f}}$) are taken as linear combinations of $\phi^{\mathrm{f}}$

$\displaystyle\rho^{\mathrm{f}}$
$\displaystyle=\rho_{l}\boldsymbol{\phi}^{\mathrm{f}}+\rho_{v}\left(1-\boldsymbol{\phi}^{\mathrm{f}}\right),$
(1) $\displaystyle\mu^{\mathrm{f}}$
$\displaystyle=\mu_{l}\boldsymbol{\phi}^{\mathrm{f}}+\mu_{v}\left(1-\boldsymbol{\phi}^{\mathrm{f}}\right),$
(2)

where $\rho_{l}$ and $\rho_{v}$ are the densities of the pure liquid and vapor
phases, respectively. $\mu_{l}$ and $\mu_{v}$ are the dynamic viscosities of
the liquid and the vapor phases.

####  2.1.1 Cavitation TEM

In the TEM approach, $\phi^{\mathrm{f}}$ is obtained as the solution of a
scalar transport equation, which can be written in the conservative form in
the ALE framework as:

$\displaystyle\left.\frac{\partial\boldsymbol{\phi}^{\mathrm{f}}}{\partial
t}\right|_{\boldsymbol{\chi}}+\boldsymbol{\phi}^{\mathrm{f}}\nabla\cdot\boldsymbol{u}^{\mathrm{f}}+\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\boldsymbol{\phi}^{\mathrm{f}}=\dfrac{\dot{m}}{\rho_{l}},$
$\displaystyle\mathrm{on}\
(\boldsymbol{x}^{\mathrm{f}},t)\in\Omega^{\mathrm{f}}$ (3)

where $\boldsymbol{\chi}$ is the referential coordinate system,
${\boldsymbol{u}^{\mathrm{f}}}={\boldsymbol{u}^{\mathrm{f}}}(\boldsymbol{x}^{\mathrm{f}},t)$
is the fluid velocity at each spatial point
$\boldsymbol{x}^{\mathrm{f}}\in\Omega^{\mathrm{f}}$ and
$\boldsymbol{u}^{\mathrm{m}}$ is the relative velocity of the spatial
coordinates $\boldsymbol{x}^{\mathrm{f}}$ with respect to the referential
coordinate system $\boldsymbol{\chi}$. The source term $\dot{m}$ is
representative of a finite mass transfer rate that governs the rates of
destruction and production of liquid by the process of cavitation. Cavitation
TEMs vary in the way $\dot{m}$ is modeled.

#### Cavitation Model A

For model A by [47], $\dot{m}$ is given as a non-linear function of
$\boldsymbol{\phi}^{\mathrm{f}}$ and $p^{\mathrm{f}}$

$\displaystyle\dot{m}_{A}=\frac{3\rho_{l}\rho_{v}}{\rho^{\mathrm{f}}R_{B}}\sqrt{\frac{2}{3\rho_{l}\left|p^{\mathrm{f}}-p_{v}\right|}}\bigg{[}C_{c}\boldsymbol{\phi}^{\mathrm{f}}$
$\displaystyle(1-\boldsymbol{\phi}^{\mathrm{f}})\operatorname{max}\left(p^{\mathrm{f}}-p_{v},0\right)$
$\displaystyle+C_{v}\boldsymbol{\phi}^{\mathrm{f}}(1+\phi_{nuc}-\boldsymbol{\phi}^{\mathrm{f}})\operatorname{min}\left(p^{\mathrm{f}}-p_{v},0\right)\bigg{]}$
(4)

Model A attempts to relate the finite mass transfer rate to the rate of
growth/collapse of an equivalent spherical bubble under an external pressure
field, using a simplification of the Rayleigh-Plesset equation. Cavitation is
assumed to initiate from nucleation sites present in the flow by a
heterogeneous nucleation process [7]. The initial concentration of nuclei per
unit volume ($n_{0}$) with an associated nuclei diameter($d_{nuc}$) and is
assumed to be a constant. It is also assumed that only vaporous cavitation
occurs, and the effect of non-condensable gases is not considered.
$R_{B}(\boldsymbol{x}^{\mathrm{f}},t)$ in Eq. (4) is representative of the
equivalent radius of the vapor volume at the coordinates
$(\boldsymbol{x}^{\mathrm{f}},t)$, while $\phi_{nuc}$ is the phase fraction of
the initial nucleation sites in an unit volume. These are calculated as

$R_{B}=\left(\frac{3}{4\pi
n_{0}}\frac{1+\phi_{nuc}-\phi^{\mathrm{f}}}{\phi^{\mathrm{f}}}\right)^{1/3}\quad\mathrm{and}\quad\phi_{nuc}=\frac{\dfrac{\pi
n_{0}d^{3}_{nuc}}{6}}{1+\dfrac{\pi n_{0}d^{3}_{nuc}}{6}}$ (5)

The vapor phase at any spatial location is assumed to be present in the form
of a concentration of bubbles with identical radii. The model requires as
input the condensation coefficient $C_{c}$ and the evaporation coefficient
$C_{v}$, which require calibration for specific flow configurations. It has
been applied to the study of different cavitating flow configurations,
including the collapse of vaporous bubbles [19] and cavitating flow over
hydrofoils [29].

We see from Eq. (4) that $\dot{m}_{A}$ depends on the local pressure
$p^{\mathrm{f}}$ as

$\displaystyle\dot{m}_{A}\propto\dfrac{p^{\mathrm{f}}-p_{v}}{\sqrt{\left|p^{\mathrm{f}}-p_{v}\right|}}$
$\displaystyle\ p^{\mathrm{f}}(\boldsymbol{x}^{\mathrm{f}},t)\in{\rm
I\\!R}-\\{0\\}$ (6)

It is observed that $\dot{m}_{A}$ is not defined when $p^{\mathrm{f}}=p_{v}$
and can lead to numerical instability. In the current work, we take

$\left.\dot{m}_{A}\right|_{p^{\mathrm{f}}=p_{v}}=\lim_{p^{\mathrm{f}}\to
p_{v}}\frac{p^{\mathrm{f}}-p_{v}}{\sqrt{\left|p^{\mathrm{f}}-p_{v}\right|}}=0$
(7)

#### Cavitation model B

For model B, $\dot{m}$ is a linear function of both
$\boldsymbol{\phi}^{\mathrm{f}}$ and the local pressure $p^{\mathrm{f}}$

$\displaystyle\dot{m}_{B}=C_{dest}\frac{\rho_{l}^{2}}{\rho_{v}}\frac{\boldsymbol{\phi}^{\mathrm{f}}}{\frac{1}{2}\rho_{l}U_{\infty}^{2}t_{\infty}}$
$\displaystyle\operatorname{min}\left(p^{\mathrm{f}}-p_{v},0\right)$
$\displaystyle+C_{prod}\rho_{l}\frac{1-\boldsymbol{\phi}^{\mathrm{f}}}{\frac{1}{2}\rho_{l}U_{\infty}^{2}t_{\infty}}\operatorname{max}\left(p^{\mathrm{f}}-p_{v},0\right)$
(8)

where $p_{v}$ is the saturation vapor pressure, $U_{\infty}$ is the free-
stream velocity and $t_{\infty}$ is the mean flow time-scale. For hydrofoils,
$t_{\infty}$ is taken as $t_{\infty}$ = $C/U_{\infty}$, where $C$ is the chord
length. $C_{dest}$ and $C_{prod}$ are semi-empirical coefficients influencing
the rates of destruction and production of liquid by the process of
cavitation. In the current work, we assume $C_{dest}=1$ and $C_{prod}=80$ [51,
25] . Model B has been derived using dimensional arguments for clusters of
large bubbles and has been applied to macro-scale cavitation over hydrofoils
[25, 6].

####  2.1.2 Fluid momentum and mass conservation

The unsteady Navier-Stokes equations for the fluid momentum and mass
conservation can be written in an ALE framework as

$\displaystyle\left.\rho^{\mathrm{f}}\frac{\partial\boldsymbol{u}^{\mathrm{f}}}{\partial
t}\right|_{\boldsymbol{\chi}}+\rho^{\mathrm{f}}\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\boldsymbol{u}^{\mathrm{f}}-\nabla\cdot\boldsymbol{\sigma}=\boldsymbol{f}^{\mathrm{f}},$
$\displaystyle\mathrm{on}\
(\boldsymbol{x}^{\mathrm{f}},t)\in\Omega^{\mathrm{f}},$ (9)
$\displaystyle\left.\frac{\partial\rho^{\mathrm{f}}}{\partial
t}\right|_{\boldsymbol{\chi}}+\rho^{\mathrm{f}}\nabla\cdot\boldsymbol{u}^{\mathrm{f}}+\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\rho^{\mathrm{f}}=0,$
$\displaystyle\mathrm{on}\
(\boldsymbol{x}^{\mathrm{f}},t)\in\Omega^{\mathrm{f}},$ (10)

where $\boldsymbol{f}^{\mathrm{f}}$ is the body force applied on the fluid and

$\boldsymbol{\sigma}=\boldsymbol{\sigma}^{\mathrm{f}}+\boldsymbol{\sigma}^{\mathrm{des}}$
(11)

where ${\boldsymbol{\sigma}^{\mathrm{f}}}$ and
$\boldsymbol{\sigma}^{\mathrm{des}}$ are the Cauchy stress tensor for a
Newtonian fluid and the turbulent stress tensor respectively, given by

$\displaystyle{\boldsymbol{\sigma}^{\mathrm{f}}}$
$\displaystyle=-{p^{\mathrm{f}}}\boldsymbol{I}+\mu^{\mathrm{f}}(\nabla{\boldsymbol{u}^{\mathrm{f}}}+(\nabla{\boldsymbol{u}^{\mathrm{f}}})^{T}),$
(12) $\displaystyle{\boldsymbol{\sigma}^{\mathrm{des}}}$
$\displaystyle=\mu_{T}(\nabla{\boldsymbol{u}^{\mathrm{f}}}+(\nabla{\boldsymbol{u}^{\mathrm{f}}})^{T}),$
(13)

where ${p^{\mathrm{f}}}$ denotes the fluid pressure and $\mu_{T}$ is the
turbulent viscosity. $\boldsymbol{\sigma}^{\mathrm{des}}$ is modeled using the
Boussinesq approximation and in the current work a hybrid URANS-LES turbulence
model is applied. The details of the turbulence model implementation can be
found in Joshi and Jaiman [32].

####  2.1.3 Convective form of TEM and local fluid compressibility

In the present work, the conservative form of the transport equation is re-
arranged in the form a convection-reaction equation. Taking the material
derivative of Eq. (1) in the ALE framework, we obtain

$\left.\frac{\partial\rho^{\mathrm{f}}}{\partial
t}\right|_{\boldsymbol{\chi}}+\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\rho^{\mathrm{f}}=\left(\rho_{l}-\rho_{v}\right)\left(\left.\frac{\partial\boldsymbol{\phi}^{\mathrm{f}}}{\partial
t}\right|_{\boldsymbol{\chi}}+\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\boldsymbol{\phi}^{\mathrm{f}}\right)$
(14)

Combining equations (3), (10) and (14), the following forms of the mass
continuity equation and the phase indicator transport equation are obtained,
which are used in the current implementation.

$\displaystyle\nabla\cdot\boldsymbol{u}^{\mathrm{f}}=\left(\frac{1}{\rho_{l}}-\frac{1}{\rho_{v}}\right)\dot{m},$
$\displaystyle\mathrm{on}\
(\boldsymbol{x}^{\mathrm{f}},t)\in\Omega^{\mathrm{f}},$ (15)
$\displaystyle\left.\frac{\partial\boldsymbol{\phi}^{\mathrm{f}}}{\partial
t}\right|_{\boldsymbol{\chi}}+\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\boldsymbol{\phi}^{\mathrm{f}}=\frac{\rho^{\mathrm{f}}}{\rho_{l}\rho_{v}}\dot{m},$
$\displaystyle\mathrm{on}\
(\boldsymbol{x}^{\mathrm{f}},t)\in\Omega^{\mathrm{f}},$ (16)

It is observed that the divergence of the velocity field is no longer zero,
and local dilation effects are introduced that are governed by the finite mass
transfer rate. This local compressibility exists only within the two-phase
mixture and the pure phases are incompressible, since for the cavitation
models chosen no mass transfer occurs when $\boldsymbol{\phi}^{\mathrm{f}}$
equals $0$ or $1$.

####  2.1.4 Turbulence viscosity modification

The URANS implementations were developed for incompressible turbulent flows,
and are not fully capable of accounting for the local compressibility induced
by the presence of vapor. For cavitating flows over hydrofoils, this leads to
an over-prediction of the turbulent stresses at the cavity closure and has
been reported by [15, 25]. This reduces the momentum of the characteristic re-
entrant jet in cavitating flows and does not allow it to penetrate the cavity.
Thus, the periodic cavity breakup and shedding is undermined and the cavity
remains attached to the surface in a quasi-steady state. This has also been
observed by the authors in preliminary investigations, although not presented
here. Hence, in the current work the turbulent viscosity is modified in the
presence of vapor. An approach similar to that adopted in [15] is taken. The
turbulent viscosity $\mu_{T}$ is modified as

$\displaystyle\mu_{T_{mod}}$
$\displaystyle=\frac{f(\rho^{\mathrm{f}})}{\rho^{\mathrm{f}}}\mu_{T}$ (17)
$\displaystyle f(\rho^{\mathrm{f}})$
$\displaystyle=\rho_{v}+\left(\frac{\rho_{v}-\rho^{\mathrm{f}}}{\rho_{v}-\rho_{l}}\right)^{n}\left(\rho_{l}-\rho_{v}\right),\quad
n\gg 1$ (18)

A value of $n=10$ is used in the current work, which has been shown to agree
well with experimental observations of cavity shedding [15].

####  2.1.5 FSI boundary conditions and fluid mesh deformation

In Section 5, we study cavitating flow over a pitching hydrofoil. We briefly
review the FSI boundary conditions and the ALE mesh motion in the continuum
setting. The modeling of FSI requires the satisfaction of the velocity
continuity and traction equilibrium at the fluid-structure boundary
$\Gamma^{\mathrm{fs}}$. Let us consider a structural domain
$\Omega^{s}\subset\mathbb{R}^{d}$ with an associated structural boundary
$\Gamma^{\mathrm{s}}(0)$ at time $t=0$. Let the mapping
$\varphi^{\mathrm{s}}\left(x^{\mathrm{s}},t\right)$ map the deformation of the
structure from its initial configuration $\Omega^{s}$ to a deformed
configuration $\Omega^{s}(t)$ at time $t$, where $x^{\mathrm{s}}$ denote the
material coordinates. We denote the initial fluid-structure interface at $t=0$
by $\Gamma^{\mathrm{fs}}(0)=\Gamma^{\mathrm{f}}(0)\cap\Gamma^{\mathrm{s}}(0)$.
At time $t$ the interface will then be deformed as
$\Gamma^{\mathrm{fs}}(t)=\varphi^{\mathrm{s}}\left(\Gamma^{\mathrm{fs}},t\right)$.
The following kinematic and dynamic conditions are satisfied on
$\Gamma^{\mathrm{fs}}$

$\displaystyle\boldsymbol{u}^{\mathrm{f}}\left(\boldsymbol{\varphi}^{\mathrm{s}}\left(\boldsymbol{x}^{\mathrm{s}},t\right),t\right)$
$\displaystyle=\boldsymbol{u}^{\mathrm{s}}\left(\boldsymbol{x}^{\mathrm{s}},t\right),$
$\displaystyle\forall\boldsymbol{x}^{\mathrm{s}}\in\Gamma^{\mathrm{fs}}$ (19)
$\displaystyle\int_{\boldsymbol{\varphi}^{\mathrm{s}}(\gamma,t)}\boldsymbol{\sigma}^{\mathrm{f}}\cdot\mathbf{n}^{\mathrm{f}}d\Gamma+\int_{\gamma}\boldsymbol{\sigma}^{\mathrm{s}}\cdot\mathbf{n}^{\mathrm{s}}d\Gamma$
$\displaystyle=0,$ $\displaystyle\forall\gamma\subset\Gamma^{\mathrm{fs}}$
(20)

where $u^{s}$ is the velocity of the structural domain,
$\mathbf{n}^{\mathrm{f}}$ and $\mathbf{n}^{\mathrm{s}}$ are the unit normals
to the deformed fluid elements $\boldsymbol{\varphi}^{\mathrm{s}}(\gamma,t)$
and their corresponding structural elements $\gamma$ on the interface
$\Gamma^{\mathrm{fs}}$ respectively. The structural stress tensor
$\sigma^{\mathrm{s}}$ is modeled depending on the type of material.

The fluid spatial coordinates are updated to conform to the structural
deformation. The motion of the coordinates which are not at
$\Gamma^{\mathrm{fs}}$ is modeled as an elastic material in equilibrium and
the mesh equation is solved as

$\displaystyle\nabla\cdot\boldsymbol{\sigma}^{\mathrm{m}}$
$\displaystyle=\mathbf{0},$ $\displaystyle\text{ on }\Omega^{\mathrm{f}},$
(21) $\displaystyle\boldsymbol{\eta}^{\mathrm{f}}$
$\displaystyle=\boldsymbol{\eta}_{D}^{\mathrm{f}},$
$\displaystyle\forall\boldsymbol{x}^{\mathrm{f}}\in\Gamma_{D}^{\mathrm{m}}$
(22)

where
$\boldsymbol{\sigma}^{\mathrm{m}}=\left(1+k_{m}\right)\left[\nabla\boldsymbol{\eta}^{\mathrm{f}}+\left(\nabla\boldsymbol{\eta}^{\mathrm{f}}\right)^{T}+\left(\nabla\cdot\boldsymbol{\eta}^{\mathrm{f}}\right)\boldsymbol{I}\right]$
is the stress experienced at the fluid spatial coordinates due to the strain
induced by the deformation of the interface, $\boldsymbol{\eta}^{\mathrm{f}}$
is the displacement of the fluid spatial coordinates. The amount of
deformation of the spatial coordinates is controlled using the local stiffness
parameter $k_{m}$. Dirichlet conditions for the fluid mesh displacement
$\boldsymbol{\eta}_{D}^{\mathrm{f}}$ are satisfied on the boundary
$\Gamma_{D}^{\mathrm{m}}$.

###  2.2 Temporal Discretization

Both the fluid mass and momentum equations and the phase indicator transport
equation are discretized in time using a generalized-$\alpha$ predictor-
corrector time integration method [14, 28]. For linear problems, the
generalized-$\alpha$ method can be second-order accurate and unconditionally
stable. This also enables the use of a single parameter called the spectral
radius $\rho_{\infty}$, to achieve user-controlled high-frequency damping.

####  2.2.1 Cavitation TEM

Let $\partial_{t}{\boldsymbol{\phi}}^{{\mathrm{f}},\mathrm{n+\alpha_{m}}}$ be
the temporal derivative of $\boldsymbol{\phi}^{\mathrm{f}}$ at time
$t^{\mathrm{n+\alpha_{m}}}$. Using the generalized-$\alpha$ method,
$\phi^{\mathrm{f}}$ is solved for at the $\mathrm{n}+\alpha$ time-level and
integrated in time $t\in\left[t^{\mathrm{n}},t^{\mathrm{n}+1}\right]$
according to the rules

$\displaystyle\phi^{\mathrm{f},\mathrm{n+1}}$
$\displaystyle=\phi^{\mathrm{f},\mathrm{n}}+\Delta
t\partial_{t}{\phi}^{\mathrm{f},\mathrm{n}}+\gamma\Delta
t(\partial_{t}{\phi}^{\mathrm{f},\mathrm{n+1}}-\partial_{t}{\phi}^{\mathrm{f},\mathrm{n}}),$
$\displaystyle\partial_{t}{\phi}^{\mathrm{f},\mathrm{n}+\alpha_{\mathrm{m}}}$
$\displaystyle=\partial_{t}{\phi}^{\mathrm{f},\mathrm{n}}+\alpha_{\mathrm{m}}(\partial_{t}{\phi}^{\mathrm{f},\mathrm{n+1}}-\partial_{t}{\phi}^{\mathrm{f},\mathrm{n}}),$
(23) $\displaystyle\phi^{\mathrm{f},\mathrm{n}+\alpha}$
$\displaystyle=\phi^{\mathrm{f},\mathrm{n}}+\alpha(\phi^{\mathrm{f},\mathrm{n+1}}-\phi^{\mathrm{f},\mathrm{n}}),$

where $\Delta t$ is the time step size and
$\partial_{t}{\phi}^{\mathrm{f},\mathrm{n}+\alpha_{\mathrm{m}}}$ is the
temporal derivative of $\phi^{\mathrm{f}}$ at the
$\mathrm{n}+\alpha_{\mathrm{m}}$ time level. $\alpha_{\mathrm{m}}$, $\alpha$
and $\gamma$ are generalized-$\alpha$ parameters based on $\rho_{\infty}$
[14].

The transport equation is arranged in the following convection-reaction form
for numerical solving

$\displaystyle
G(\partial_{t}{\phi}^{\mathrm{f},\mathrm{n}+\alpha_{\mathrm{m}}},\phi^{\mathrm{f},\mathrm{n}+\alpha})=\partial_{t}{\phi}^{\mathrm{f},\mathrm{n}+\alpha_{\mathrm{m}}}+\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\phi^{\mathrm{f},\mathrm{n}+\alpha}+s\phi^{\mathrm{f},\mathrm{n}+\alpha}-f=0,$
(24)

where $\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}$ is the
convection velocity. The reaction coefficient $s$ and the source term $f$ for
model A are given by

$\displaystyle s_{A}$
$\displaystyle=-\frac{3}{R_{B}}\sqrt{\frac{2}{3\rho_{l}\left|p^{\mathrm{f}}-p_{v}\right|}}\bigg{[}C_{c}(1-\boldsymbol{\phi}^{\mathrm{f}})\operatorname{max}\left(p^{\mathrm{f}}-p_{v},0\right)+C_{v}(1+\phi_{nuc}-\boldsymbol{\phi}^{\mathrm{f}})\operatorname{min}\left(p^{\mathrm{f}}-p_{v},0\right)\bigg{]},$
(25) $\displaystyle f_{A}$ $\displaystyle=0.$ (26)

while for model B, we have

$\displaystyle s_{B}$
$\displaystyle=\frac{\rho^{\mathrm{f}}}{\rho_{v}}\bigg{[}-\frac{\rho_{l}}{\rho_{v}}\frac{C_{dest}}{\frac{1}{2}\rho_{l}U_{\infty}^{2}t_{\infty}}\operatorname{min}\left(p-p_{v},0\right)+\frac{C_{prod}}{\frac{1}{2}\rho_{l}U_{\infty}^{2}t_{\infty}}\operatorname{max}\left(p-p_{v},0\right)\bigg{]},$
(27) $\displaystyle f_{B}$
$\displaystyle=\frac{\rho^{\mathrm{f}}}{\rho_{v}}\bigg{[}\frac{C_{prod}}{\frac{1}{2}\rho_{l}U_{\infty}^{2}t_{\infty}}\operatorname{max}\left(p-p_{v},0\right)\bigg{]}.$
(28)

####  2.2.2 Navier-Stokes equations

The Navier-Stokes equations are also discretized using the
generalized-$\alpha$ time integration for consistency with the phase indicator
transport equation. The variational formulation employs the following
relations for time integration during
$t\in\left[t^{\mathrm{n}},t^{\mathrm{n}+1}\right]$

$\displaystyle\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}$
$\displaystyle=\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\Delta
t\partial_{t}\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\gamma\Delta
t(\partial_{t}\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}-\partial_{t}\boldsymbol{u}^{\mathrm{f},\mathrm{n}}),$
(29) $\displaystyle\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}$
$\displaystyle=\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\alpha(\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}-\boldsymbol{u}^{\mathrm{f},\mathrm{n}}),$
(30)
$\displaystyle\partial_{t}\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha_{m}}}$
$\displaystyle=\partial_{t}\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\alpha_{\mathrm{m}}(\partial_{t}\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}-\partial_{t}\boldsymbol{u}^{\mathrm{f},\mathrm{n}}).$
(31)

The semi-discrete forms of the Navier-Stokes equations are written as

$\displaystyle\left.\rho^{\mathrm{f}}\partial_{t}\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha_{\mathrm{m}}}\right|_{\chi}+\rho^{\mathrm{f}}\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m},\mathrm{n}+\alpha}\right)\cdot\nabla\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\nabla\cdot\boldsymbol{\sigma}^{\mathrm{n}+\alpha}-\boldsymbol{f}^{\mathrm{n}+\alpha}$
$\displaystyle=0,$ (32)
$\displaystyle\nabla\cdot\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\left(\frac{1}{\rho_{l}}-\frac{1}{\rho_{v}}\right)\dot{m}$
$\displaystyle=0,$ (33)

###  2.3 Spatial discretization and variational statement

Next, we present the stabilized variational statements of the governing
equations. We first present the positivity preserving variational form of the
cavitation model, followed by the Navier-Stokes equations. The stabilized
variational form for the incompressible Navier-Stokes has been discussed in
detail in other work [32, 31], which we shall not reproduce here. However, the
presence of the non-zero divergence of the fluid velocity introduces
additional terms to the formulation, which merits a brief discussion.

####  2.3.1 Cavitation TEM

The fluid computational domain $\Omega^{\mathrm{f}}$ is spatially discretized
into $\mathrm{n_{el}}$ number of elements such that
$\Omega^{\mathrm{f}}=\cup_{\mathrm{e=1}}^{\mathrm{n_{el}}}\Omega^{{\mathrm{f}},\mathrm{e}}$
and
$\emptyset=\cap_{\mathrm{e=1}}^{\mathrm{n_{el}}}\Omega^{{\mathrm{f}},\mathrm{e}}$.
The trial solutions are taken from from the space $\mathcal{S}^{\mathrm{h}}$,
which equal the given Dirichlet boundary condition at the boundary
$\Gamma_{D}$. The test functions are taken from the space
$\mathcal{V}^{\mathrm{h}}$, which vanish on the Dirichlet boundary. We state
the variational form of the phase indicator transport equation to find
$\phi^{\mathrm{f}}_{\mathrm{h}}(\boldsymbol{x}^{\mathrm{f}},t^{\mathrm{n}+\alpha})\in\mathcal{S}^{\mathrm{h}}$
such that $\forall w_{\mathrm{h}}\in\mathcal{V}^{\mathrm{h}}$,

$\displaystyle\int_{\Omega^{\mathrm{f}}}\bigg{(}w_{\mathrm{h}}\partial_{t}{\phi}^{\mathrm{f},\mathrm{n+\alpha_{m}}}_{\mathrm{h}}+w_{\mathrm{h}}\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m},\mathrm{n}+\alpha}\right)\cdot\nabla\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}+w_{\mathrm{h}}s\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}-w_{\mathrm{h}}f\bigg{)}\mathrm{d}\Omega^{\mathrm{f}}$
$\displaystyle+\displaystyle\sum_{\mathrm{e=1}}^{\mathrm{n_{el}}}\int_{\Omega^{{\mathrm{f}},\mathrm{e}}}\bigg{(}\big{(}\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla
w_{\mathrm{h}}\big{)}\tau_{\phi}\big{(}\partial_{t}{\phi^{\mathrm{f},\mathrm{n}+\alpha_{m}}}_{\mathrm{h}}+\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}+s\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}-f\big{)}\bigg{)}\mathrm{d}\Omega^{{\mathrm{f}},\mathrm{e}}$
$\displaystyle+\displaystyle\sum_{\mathrm{e=1}}^{\mathrm{n_{el}}}\int_{\Omega^{{\mathrm{f}},\mathrm{e}}}\chi\frac{|\boldsymbol{\mathcal{R}}_{\phi}|}{|\nabla\phi^{\mathrm{f}}_{\mathrm{h}}|}k_{s}^{\mathrm{add}}\nabla
w_{\mathrm{h}}\cdot\bigg{(}\frac{\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)\otimes\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)}{|\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)|^{2}}\bigg{)}\cdot\nabla\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}\mathrm{d}\Omega^{{\mathrm{f}},\mathrm{e}}$
(34)
$\displaystyle+\sum_{\mathrm{e=1}}^{\mathrm{n_{el}}}\int_{\Omega^{{\mathrm{f}},\mathrm{e}}}\chi\frac{|\boldsymbol{\mathcal{R}}_{\phi}|}{|\nabla\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}|}k^{\mathrm{add}}_{c}\nabla
w_{\mathrm{h}}\cdot\bigg{(}\mathbf{I}-\frac{\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)\otimes\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)}{|\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)|^{2}}\bigg{)}\cdot\nabla\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}\mathrm{d}\Omega^{{\mathrm{f}},\mathrm{e}}$
$\displaystyle=0,$

where the first line represents the standard Galerkin finite element terms.
The second line consists of linear stabilization terms with the stabilization
parameter $\tau_{\phi}$ given by [53]

$\displaystyle\tau_{\phi}$ $\displaystyle=\bigg{[}\bigg{(}\frac{2}{\Delta
t}\bigg{)}^{2}+\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\boldsymbol{G}\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)+s^{2}\bigg{]}^{-1/2},$
(35)

where $\boldsymbol{G}$ is the element contravariant metric tensor, which is
defined as

$\displaystyle\boldsymbol{G}=\frac{\partial\boldsymbol{\xi}^{T}}{\partial\boldsymbol{x}^{\mathrm{f}}}\frac{\partial\boldsymbol{\xi}}{\partial\boldsymbol{x}^{\mathrm{f}}},$
(36)

where $\boldsymbol{x}^{\mathrm{f}}$ and $\boldsymbol{\xi}$ are the physical
and parametric coordinates respectively. $\boldsymbol{\mathcal{R}}_{\phi}$ is
the residual of the phase indicator transport equation given as

$\displaystyle\boldsymbol{\mathcal{R}}_{\phi}=\partial_{t}\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}+\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}+s\phi^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}-f.$
(37)

The linear stabilization terms are used to address spurious oscillations in
the solution when convection and reaction effects are dominant. However,
separate treatment is required to address oscillations in the solution near
the regions of high gradients [33]. Across the cavity interface, the
transition from the liquid to the vapor phase takes place over the span of a
few elements, and the spatial gradient of $\phi^{\mathrm{f}}$ is very high. As
seen in Eq. (1) and Eq. (2), physical properties of the fluid such as the
density and the viscosity are obtained as weighted linear interpolations of
$\phi^{\mathrm{f}}$. Unbounded oscillations in $\phi^{\mathrm{f}}$ can result
in negative values of $\rho^{\mathrm{f}}$ and $\mu^{\mathrm{f}}$, which are
unphysical and can induce numerical instability. To address this, we introduce
additional non-linear stabilization terms that impart a positivity property to
the underlying element-level matrices. The third and fourth lines of Eq.
(2.3.1) contain the positivity preserving nonlinear stabilization terms in the
streamwise and crosswind directions respectively[33]. These essentially act as
added diffusion in the region of high gradients and ensure that the element-
level matrix is an $M$-matrix.

The PPV parameters $\chi$, $k_{s}^{\mathrm{add}}$ and $k_{c}^{\mathrm{add}}$
for the phase indicator transport equation are obtained as:

$\displaystyle\chi$
$\displaystyle=\frac{2}{|s|h+2|\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)|},$
(38) $\displaystyle k_{s}^{\mathrm{add}}$
$\displaystyle=\mathrm{max}\bigg{\\{}\frac{||\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)|-\tau_{\phi}|\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)|s|h}{2}-\tau_{\phi}|\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)|^{2}+\frac{sh^{2}}{6},0\bigg{\\}},$
(39) $\displaystyle k_{c}^{\mathrm{add}}$
$\displaystyle=\mathrm{max}\bigg{\\{}\frac{|\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)|h}{2}+\frac{sh^{2}}{6},0\bigg{\\}},$
(40)

where $|\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)|$
is the magnitude of the convection velocity and $h$ is the characteristic
element length [31]. We next present the Navier-Stokes equations and its
variational form for the two-phase solver.

####  2.3.2 Navier-Stokes equations

The trial solution is taken from the function space
$\mathcal{S}^{\mathrm{h}}$, the values of which satisfy the Dirichlet boundary
condition at the Dirichlet boundary $\Gamma_{D}$. Let
$\mathcal{V}^{\mathrm{h}}$ be the space of test functions which vanish on
$\Gamma_{D}$. We formulate the variational statement for the fluid flow
equations to find
$[{\boldsymbol{u}}_{\mathrm{h}}^{\mathrm{n+\alpha}},{p}_{\mathrm{h}}^{\mathrm{n+1}}]\in\mathcal{S}^{\mathrm{h}}$
such that
$\forall[\boldsymbol{\psi}_{\mathrm{h}},q_{\mathrm{h}}]\in\mathcal{V}^{\mathrm{h}}$,

$\displaystyle\int_{\Omega^{\mathrm{f}}}\boldsymbol{\psi}_{h}^{\mathrm{f}}\cdot\left(\left.\rho^{\mathrm{f}}\partial_{t}\boldsymbol{u}_{h}^{\mathrm{f},\mathrm{n}+\alpha_{\mathrm{m}}}\right|_{\boldsymbol{\chi}}+\rho^{\mathrm{f}}\left(\boldsymbol{u}_{h}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\boldsymbol{u}_{h}^{\mathrm{f},\mathrm{n}+\alpha}\right)d\Omega$
$\displaystyle+\int_{\Omega^{\mathrm{f}}}\nabla\boldsymbol{\psi}_{h}^{\mathrm{f}}:\boldsymbol{\sigma}_{h}^{\mathrm{n}+\alpha}d\Omega+\int_{\Omega^{\mathrm{f}}\left(t^{\mathrm{n}+1}\right)}q_{h}\left(\nabla\cdot\boldsymbol{u}_{h}^{\mathrm{f},\mathrm{n}+\alpha}\right)d\Omega$
$\displaystyle+\sum_{e=1}^{\mathrm{n}_{\mathrm{el}}}\int_{\Omega^{\mathrm{f}}}\frac{\tau_{m}}{\rho^{\mathrm{f}}}\left(\rho^{\mathrm{f}}\left(\boldsymbol{u}_{h}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\boldsymbol{\psi}_{h}^{\mathrm{f}}+\nabla
q_{h}\right)\cdot\boldsymbol{\mathcal{R}}_{m}d\Omega^{\mathrm{e}}+\sum_{e=1}^{\mathrm{n}_{\mathrm{el}}}\int_{\Omega^{e}}\nabla\cdot\boldsymbol{\psi}_{h}^{\mathrm{f}}\tau_{c}\rho^{\mathrm{f}}\boldsymbol{\mathcal{R}}_{c}d\Omega^{e}$
(41)
$\displaystyle-\displaystyle\sum_{\mathrm{e=1}}^{\mathrm{n_{el}}}\int_{\Omega^{\mathrm{e}}}\tau_{\mathrm{m}}\boldsymbol{\psi}_{\mathrm{h}}\cdot(\boldsymbol{\mathcal{R}}_{\mathrm{m}}\cdot\nabla{\boldsymbol{u}}_{\mathrm{h}}^{\mathrm{n+\alpha}})\mathrm{d\Omega^{e}}-\displaystyle\sum_{\mathrm{e=1}}^{\mathrm{n_{el}}}\int_{\Omega^{\mathrm{e}}}\frac{\nabla\boldsymbol{\psi}_{\mathrm{h}}}{\rho(\phi)}:(\tau_{\mathrm{m}}\boldsymbol{\mathcal{R}}_{\mathrm{m}}\otimes\tau_{\mathrm{m}}\boldsymbol{\mathcal{R}}_{\mathrm{m}})\mathrm{d\Omega^{e}}$
$\displaystyle=\int_{\Omega^{\mathrm{f}}}\boldsymbol{\psi}_{h}^{\mathrm{f}}\cdot\boldsymbol{f}^{\mathrm{f},\mathrm{n}+\alpha}\;d\Omega+\int_{\Gamma_{N}^{\mathrm{f}}}\boldsymbol{\psi}_{h}^{\mathrm{f}}\cdot\boldsymbol{h}^{\mathrm{f},\mathrm{n}+\alpha}\;d\Gamma+\underbrace{\int_{\Omega^{\mathrm{f}}}q_{h}\left(\frac{1}{\rho_{l}}-\frac{1}{\rho_{v}}\right)\dot{m}^{\mathrm{f},\mathrm{n}+\alpha}\;d\Omega}$

where the first and second lines contain the standard Galerkin finite element
terms of the momentum and the continuity equation. The third line contains the
Galerkin Least Squares stabilization terms for the momentum and mass
continuity equations. The fourth line contains stabilization terms based on
the multi-scale argument [26, 24]. The fifth line contains the Galerkin terms
for the body force and the Neumann boundary in the momentum equation. The term
in under-braces corresponds to the Galerkin projection of the term in the mass
conservation equation dependent on $p^{\mathrm{f}}$, and has been introduced
in this formulation. The element-wise residual of the momentum and the
continuity equations denoted by $\boldsymbol{\mathcal{R}}_{\mathrm{m}}$ and
$\boldsymbol{\mathcal{R}}_{\mathrm{c}}$ respectively are given by

$\displaystyle\boldsymbol{\mathcal{R}}_{\mathrm{m}}({\boldsymbol{u}^{\mathrm{f}}},{p}^{\mathrm{f}})$
$\displaystyle=\rho^{\mathrm{f}}\partial_{t}{\boldsymbol{u}}_{\mathrm{h}}^{\mathrm{n+\alpha_{m}}}+\rho^{\mathrm{f}}\left(\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla{\boldsymbol{u}}_{\mathrm{h}}^{\mathrm{f},\mathrm{n}+\alpha}-\nabla\cdot{\boldsymbol{\sigma}}_{\mathrm{h}}^{\mathrm{f},\mathrm{n}+\alpha}-\boldsymbol{f}^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}},$
(42)
$\displaystyle\boldsymbol{\mathcal{R}}_{\mathrm{c}}({\boldsymbol{u}^{\mathrm{f}}},{p}^{\mathrm{f}},\phi^{\mathrm{f}})$
$\displaystyle=\nabla\cdot\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}_{\mathrm{h}}-\left(\frac{1}{\rho_{l}}-\frac{1}{\rho_{v}}\right)\dot{m}^{\mathrm{f},\mathrm{n}+\alpha}.$
(43)

$\tau_{\mathrm{m}}$ and $\tau_{\mathrm{c}}$ are stabilization parameters [8,
57, 18] defined as

$\displaystyle\tau_{\mathrm{m}}$
$\displaystyle=\bigg{[}\bigg{(}\frac{2}{\Delta
t}\bigg{)}^{2}+{\boldsymbol{u}}_{\mathrm{h}}\cdot\boldsymbol{G}{\boldsymbol{u}}_{\mathrm{h}}+C_{I}\bigg{(}\frac{\mu(\phi)}{\rho(\phi)}\bigg{)}^{2}\boldsymbol{G}:\boldsymbol{G}\bigg{]}^{-1/2},$
(44) $\displaystyle\tau_{\mathrm{c}}$
$\displaystyle=\frac{1}{\mathrm{tr}(\boldsymbol{G})\tau_{\mathrm{m}}},$ (45)

where $C_{I}$ is a constant derived from the element-wise inverse estimate
[23] and $\mathrm{tr(\boldsymbol{G})}$ denotes the trace of $\boldsymbol{G}$.

###  2.4 Finite element matrix form

The Navier-Stokes equations are coupled with the phase fraction transport
equation in a partitioned manner as shown in Eq. (46) and Eq. (47). The
equations are linearized with a Newton-Raphson technique to solve for the
incremental velocity, pressure and phase indicator which are advanced in time
using the Generalized-$\alpha$ time integration method. The linearized
matrices for the fluid flow and the phase indicator transport equations are
arranged in the form

$\left[\begin{array}[]{cc}\boldsymbol{K}_{\Omega^{\mathrm{f}}}&\boldsymbol{G}_{\Omega^{\mathrm{f}}}\\\
-\boldsymbol{G}_{\Omega^{\mathrm{f}}}^{T}&\boldsymbol{C}_{\Omega^{\mathrm{f}}}\end{array}\right]\left\\{\begin{array}[]{c}\Delta\boldsymbol{u}^{\mathrm{f},n+\alpha}\\\
\Delta
p^{\mathrm{f},n+1}\end{array}\right\\}=\left\\{\begin{array}[]{l}{\boldsymbol{\mathcal{R}}}_{\mathrm{m}}\\\
{\boldsymbol{\mathcal{R}}}_{\mathrm{c}}\end{array}\right\\}$ (46)
$\left[\boldsymbol{K}_{\phi}\right]\\{\Delta\phi^{\mathrm{f},n+\alpha}\\}=\\{{\boldsymbol{\mathcal{R}}}_{\phi}\\}$
(47)

where $\boldsymbol{K}_{\Omega^{\mathrm{f}}}$ is the stiffness matrix of the
momentum equation, $\boldsymbol{G}_{\Omega^{\mathrm{f}}}$ is the discrete
gradient operator and $\boldsymbol{G}^{T}_{\Omega^{\mathrm{f}}}$ is the
divergence operator. $\boldsymbol{C}_{\Omega^{\mathrm{f}}}$ consists of the
pressure-pressure stabilization term and the terms in the mass continuity
equation having dependency on the pressure. $\boldsymbol{K}_{\phi}$ is the
stiffness matrix for the phase indicator transport equation, consisting of the
transient, convection, reaction, linear stabilization terms and the non-linear
PPV stabilization terms. $\Delta\boldsymbol{u}^{\mathrm{f},n+\alpha}$, $\Delta
p^{\mathrm{f},n+1}$ and $\Delta\phi^{\mathrm{f},n+\alpha}$ are the increments
in velocity, pressure and the phase indicator respectively.
Linearization of the phase indicator transport equation requires the
derivative of the terms in Eq. (2.2.1) with respect to
$\phi^{\mathrm{f},n+\alpha}$.

$\displaystyle\dfrac{\partial G}{\partial\phi^{\mathrm{f},n+\alpha}}$
$\displaystyle=\dfrac{\partial\left(\partial_{t}{\phi}^{\mathrm{f},\mathrm{n}+\alpha_{\mathrm{m}}}\right)}{\partial\phi^{\mathrm{f},n+\alpha}}+\cancelto{0}{\dfrac{\partial\left(\left(\boldsymbol{u}^{\mathrm{f}}-\boldsymbol{u}^{\mathrm{m}}\right)\cdot\nabla\phi^{\mathrm{f},\mathrm{n}+\alpha}\right)}{\partial\phi^{\mathrm{f},n+\alpha}}}+\dfrac{\partial\left(s\phi^{\mathrm{f},\mathrm{n}+\alpha}\right)}{\partial\phi^{\mathrm{f},n+\alpha}}-\cancelto{0}{\dfrac{\partial
f}{\partial\phi^{\mathrm{f},n+\alpha}}}$ (48)
$\displaystyle=\dfrac{\alpha_{\mathrm{m}}}{\alpha\gamma\Delta
t}+\dfrac{\partial\left(s\phi^{\mathrm{f},\mathrm{n}+\alpha}\right)}{\partial\phi^{\mathrm{f},n+\alpha}}$
(49)

where the derivative of the transient term reduces to a constant obtained from
the generalized-$\alpha$ relations in Eq. (23). Owing to the symmetrical
property of second order cross-derivatives, the convective term also reduces
to zero. The derivative of the reaction term
$\dfrac{\partial\left(s\phi^{\mathrm{f},\mathrm{n}+\alpha}\right)}{\partial\phi^{\mathrm{f},\mathrm{n}+\alpha}}$
needs to be calculated. Similarly, in the linearization of the mass continuity
equation we encounter the term
$\dfrac{\partial\dot{m}^{\mathrm{f},\mathrm{n}+\alpha}}{\partial
p^{\mathrm{f},\mathrm{n}+1}}$. Care must be taken in forming the terms
$\dfrac{\partial\dot{m}}{\partial p^{\mathrm{f},\mathrm{n}+1}}$ and
$\dfrac{\partial\left(s\phi^{\mathrm{f},\mathrm{n}+\alpha}\right)}{\partial\phi^{\mathrm{f},\mathrm{n}+\alpha}}$
for a stable linearization of the matrices
$\boldsymbol{C}_{\Omega^{\mathrm{f}}}$ and $\boldsymbol{K}_{\phi}$
respectively. In model A, $\dot{m}$ varies non-linearly with both $p$ and
$\phi^{\mathrm{f}}$, along with a discontinuity at $p=p_{v}$. Model B presents
a simpler linearization with $\dot{m}$ being a linear function of both $p$ and
$\phi^{\mathrm{f}}$. We propose linearlizations for the said terms in Tables.
(1) and (2) respectively. The presented linearlizations have been tested to be
stable over a range of temporal and spatial discretizations on two different
cavitating flow configurations.

$\Omega^{\mathrm{f}}$$\Omega^{\mathrm{s}}$$\Gamma^{\mathrm{fs}}$collapsing cavitiesresonating tip vortex cavityre-entrant jetshed cavitiesleading-edge partial cavities$\&$ vortices$U_{\infty}$ |  Model A |  Model B
---|---|---
$p^{\mathrm{n}+1}_{h}>p_{v}$ |  $C_{c}\boldsymbol{\phi}^{\mathrm{f},\mathrm{n}+\alpha}_{h}(1-\boldsymbol{\phi}^{\mathrm{f},\mathrm{n}+\alpha}_{h})\dfrac{\rho_{l}\rho_{v}}{\rho^{\mathrm{f}}R_{B}}\sqrt{\dfrac{3}{2\rho_{l}\left(p^{\mathrm{n}+1}_{h}-p_{v}\right)}}$ |  $C_{prod}\rho_{l}\dfrac{1-\boldsymbol{\phi}^{\mathrm{f},\mathrm{n}+\alpha}_{h}}{\dfrac{1}{2}\rho_{l}U_{\infty}^{2}t_{\infty}}$
$p^{\mathrm{n}+1}_{h}<p_{v}$ |  $C_{v}\boldsymbol{\phi}^{\mathrm{f},\mathrm{n}+\alpha}_{h}(1+\phi_{nuc}-\boldsymbol{\phi}^{\mathrm{f},\mathrm{n}+\alpha}_{h})\dfrac{\rho_{l}\rho_{v}}{\rho^{\mathrm{f}}R_{B}}\sqrt{\dfrac{3}{2\rho_{l}\left(p_{v}-p^{\mathrm{n}+1}_{h}\right)}}$ |  $C_{dest}\dfrac{\rho_{l}^{2}}{\rho_{v}}\dfrac{\boldsymbol{\phi}^{\mathrm{f},\mathrm{n}+\alpha}_{h}}{\dfrac{1}{2}\rho_{l}U_{\infty}^{2}t_{\infty}}$
$p^{\mathrm{n}+1}_{h}=p_{v}$ |  0, using $\lim_{p^{\mathrm{n}+1}_{h}\to p_{v}}\dfrac{p^{\mathrm{n}+1}_{h}-p_{v}}{\sqrt{\left|p^{\mathrm{n}+1}_{h}-p_{v}\right|}}=0$ |  0

Table 1: $\dfrac{\partial\dot{m}^{\mathrm{f},\mathrm{n}+\alpha}}{\partial
p^{\mathrm{f},\mathrm{n}+1}}$ in linearized matrix
$\boldsymbol{C}_{\Omega^{\mathrm{f}}}$

$\Omega^{\mathrm{f}}$$\Omega^{\mathrm{s}}$$\Gamma^{\mathrm{fs}}$collapsing cavitiesresonating tip vortex cavityre-entrant jetshed cavitiesleading-edge partial cavities$\&$ vortices$U_{\infty}$ |  Model A |  Model B
---|---|---
$p^{\mathrm{n}+1}_{h}>p_{v}$ |  $-C_{c}(1-2\boldsymbol{\phi}^{\mathrm{f},\mathrm{n}+\alpha}_{h})\dfrac{3}{R_{B}}\sqrt{\dfrac{2\left(p^{\mathrm{n}+1}_{h}-p_{v}\right)}{3\rho_{l}}}$ |  $\dfrac{\rho^{\mathrm{f}}}{\rho_{v}}\dfrac{C_{prod}}{\dfrac{1}{2}\rho_{l}U_{\infty}^{2}t_{\infty}}\left(p^{\mathrm{n}+1}_{h}-p_{v}\right)$
$p^{\mathrm{n}+1}_{h}<p_{v}$ |  $C_{v}(1+\phi_{nuc}-2\boldsymbol{\phi}^{\mathrm{f},\mathrm{n}+\alpha}_{h})\dfrac{3}{R_{B}}\sqrt{\dfrac{2\left(p_{v}-p^{\mathrm{n}+1}_{h}\right)}{3\rho_{l}}}$ |  $-\dfrac{\rho^{\mathrm{f}}\rho_{l}}{\rho^{2}_{v}}\dfrac{C_{dest}}{\dfrac{1}{2}\rho_{l}U_{\infty}^{2}t_{\infty}}\left(p^{\mathrm{n}+1}_{h}-p_{v}\right)$
$p^{\mathrm{n}+1}_{h}=p_{v}$ |  0, using $\lim_{p^{\mathrm{n}+1}_{h}\to p_{v}}\dfrac{p^{\mathrm{n}+1}_{h}-p_{v}}{\sqrt{\left|p^{\mathrm{n}+1}_{h}-p_{v}\right|}}=0$ |  0

Table 2:
$\dfrac{\partial\left(s\phi^{\mathrm{f},\mathrm{n}+\alpha}\right)}{\partial\phi^{\mathrm{f},\mathrm{n}+\alpha}}$
in linearized matrix $\boldsymbol{K}_{\phi}$

###  2.5 Implementation details

Algorithm 1 shows the algorithm for the staggered partitioned coupling of the
implicit Navier-Stokes and cavitation solvers. This provides two unique
advantages compared to monolithic couplings - (i) simplicity of linearization
and implementation, and (ii) the flexibility to couple additional governing
equations as demanded by the case being studied. For example, in Sections 4
and 5, the additional equations for turbulence modeling and the ALE mesh
update are similarly coupled using a staggered partitioned approach. For
stable and accurate convergence, non-linear predictor-multicorrector
iterations are performed within each time step. Let us consider the time level
$t^{\mathrm{n}}$ and the associated velocity
$\boldsymbol{u}^{\mathrm{f}}(\boldsymbol{x},t^{\mathrm{n}})$, pressure
$p^{\mathrm{f}}(\boldsymbol{x},t^{\mathrm{n}})$ and phase indicator
$\phi^{\mathrm{f}}(\boldsymbol{x},t^{\mathrm{n}})$ fields. Each non-linear
predictor-corrector iteration consists of one pass through the steps [A], [B],
[C] and [D]. In step [A] of the predictor-corrector iteration $\mathrm{k}$, a
predictor fluid velocity
$\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}$ and pressure
$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}$ is obtained from the solution
of the Navier-Stokes equations(Eq. (46)). These are passed to the cavitation
solver in step [B]. An updated phase indicator
$\phi^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}$ is obtained in step [C] by
solving Eq. (47), which is then passed to the Navier-Stokes solver in step
[D]. This updated phase indicator value is used to interpolate the fluid
density and dynamic viscosity, and prepare the Navier-Stokes matrices for the
next iteration $\mathrm{k+1}$. This cyclic process is continued till the
solvers have achieved the convergence criteria. Then the coupled solver is
advanced to the next time level $t^{\mathrm{n+1}}$.

Algorithm 1 Partitioned coupling of implicit Navier-Stokes and cavitation
solvers

Input: $\boldsymbol{u}^{\mathrm{f},0}$, $p^{\mathrm{f},0}$,
$\phi^{\mathrm{f},0}$
for n $\leftarrow$ 0 to $n_{last}$ do

Predict solution
---
$\bigg{[}\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}_{(0)}\;\;p^{\mathrm{f},\mathrm{n+1}}_{(0)}\;\;\phi^{\mathrm{f},\mathrm{n+1}}_{(0)}\bigg{]}\leftarrow\bigg{[}\boldsymbol{u}^{\mathrm{f},\mathrm{n}}\;\;p^{\mathrm{f},\mathrm{n}}\;\;\phi^{\mathrm{f},\mathrm{n}}\bigg{]}$
Interpolate density and viscosity fields
$\rho^{\mathrm{f}}(\phi^{\mathrm{f},\mathrm{n}+1}_{(0)})$,
$\mu^{\mathrm{f}}(\phi^{\mathrm{f},\mathrm{n}+1}_{(0)})$
for k $\leftarrow$ 0 to convergence/$k_{max}$ do
$\Omega^{\mathrm{f}}$$\Omega^{\mathrm{s}}$$\Gamma^{\mathrm{fs}}$collapsing
cavitiesresonating tip vortex cavityre-entrant jetshed cavitiesleading-edge
partial cavities$\&$ vortices$U_{\infty}$[A] Navier-Stokes Solver[C]
Cavitation SolverNon-
linear[D]$\boldsymbol{u}^{\mathrm{n+1}}_{\mathrm{(k+1)}},p^{\mathrm{n+1}}_{\mathrm{(k+1)}}$1\.
Interpolate solution2\. Solve3\. Correct solution4\. Update solution1\.
Interpolate solution2\. Solve3\. Correct solution4\. Update
solution$\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\alpha(\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k)}}-\boldsymbol{u}^{\mathrm{f},\mathrm{n}})$$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow
p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k)}}$$\Delta\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}$
and $\Delta p^{\mathrm{f},\mathrm{n+1}}$ in Eq.
(46)$\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}+\Delta\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}$$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow
p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}+\Delta
p^{\mathrm{f},\mathrm{n+1}}$$\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\dfrac{1}{\alpha}(\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}-\boldsymbol{u}^{\mathrm{f},\mathrm{n}})$$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow
p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}+\Delta
p^{\mathrm{f},\mathrm{n+1}}$$\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\phi^{\mathrm{f},\mathrm{n}}+\alpha(\phi^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k)}}-\phi^{\mathrm{f},\mathrm{n}})$$\Delta\phi^{\mathrm{f},\mathrm{n+\alpha}}$
in Eq.
(47)$\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}+\Delta\phi^{\mathrm{f},\mathrm{n+\alpha}}$$\phi^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow\phi^{\mathrm{f},\mathrm{n}}+\dfrac{1}{\alpha}(\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}-\phi^{\mathrm{f},\mathrm{n}})$iterations$\phi^{\mathrm{n+1}}_{\mathrm{(k+1)}}$[B]
end do

end do

A similar partitioned iterative coupling is used between the ALE fluid-
cavitation and the structural displacement. In the current study, the
structural displacement is dictated by the prescribed motion. Alternatively,
for fully-coupled FSI studies the structural displacement can be obtained by
solving the governing equations for the structural mechanics as presented in
Joshi and Jaiman [32]). At time $t^{\mathrm{n}}$, let us consider a predictor
structural displacement
$\boldsymbol{\eta}^{\mathrm{s}}\left(\boldsymbol{x}^{\mathrm{s}},t^{\mathrm{n}}\right)$
obtained from the structural update in the Lagrangian reference frame. These
structural displacements are then passed to the fluid solvers, while
satisfying kinematic conditions at the fluid-structure interface
$\Gamma^{\mathrm{fs}}$ as follows. At the time $t^{\mathrm{n+1}}$, the
displacement $\boldsymbol{\eta}^{\mathrm{f}}$ of the fluid spatial coordinates
(i,e., Eulerian mesh coordinates) at the wetted boundary
$\Gamma^{\mathrm{fs}}$ are equated to the structural displacement
$\boldsymbol{\eta}^{\mathrm{s}}$.

$\displaystyle\boldsymbol{\eta}^{\mathrm{f},\mathrm{n}+1}=\boldsymbol{\eta}^{\mathrm{s}},$
$\displaystyle\text{ on }\Gamma^{\mathrm{fs}}$ (50)

In addition, the velocity continuity is satisfied at the interface is
satisfied at the time $t^{\mathrm{n}+\alpha}$ as

$\displaystyle\boldsymbol{u}^{\mathrm{f},\mathrm{n}+\alpha}=\boldsymbol{u}^{\mathrm{m},\mathrm{n}+\alpha},$
$\displaystyle\text{ on }\Gamma^{\mathrm{fs}}$ (51)

with the mesh velocity at the interface evaluated as

$\displaystyle\boldsymbol{u}^{\mathrm{m},\mathrm{n}+\alpha}=\frac{\boldsymbol{\eta}^{\mathrm{f},\mathrm{n}+1}-\boldsymbol{\eta}^{\mathrm{f},\mathrm{n}}}{\Delta
t},$ $\displaystyle\text{ on }\Gamma^{\mathrm{fs}}$ (52)

Away from the interface $\Gamma^{\mathrm{fs}}$ and any Dirichlet conditions on
$\boldsymbol{\eta}^{\mathrm{f}}$ on the Dirichlet boundary
$\Gamma_{D}^{\mathrm{m}}$, Eq. (21) is solved for updating the fluid spatial
coordinates. The flow-cavitation equations are then solved with the updated
kinematic boundary conditions and the displaced fluid spatial coordinates. The
updated hydrodynamic forces are passed to the structural solver to correct the
deformations. This sequence is repeated iteratively till converged solutions
are obtained.

For the finite element discretization of the variables
$\boldsymbol{u}^{\mathrm{f}}$, $p^{\mathrm{f}}$ and $\phi^{\mathrm{f}}$, we
consider equal-order interpolations. The Harwell-Boeing sparse matrix format
is used to form and store the matrices for the linear system of equations. A
Generalized Minimal RESidual (GMRES)[46] algorithm is used to solve the linear
system. We observe that 2-3 non-linear iterations(as described in Algorithm
(1)) are sufficient to obtain a converged solution at each time-step. The
solver uses communication protocols based on standard message passing
interface [1] for parallel computing on distributed memory clusters.

##  3 Numerical verification and convergence

A stabilized numerical method for the solution of cavitating flows has been
presented. In this section we verify the implementation and discuss its
convergence and efficacy.

###  3.1 Analytical solution of vaporous spherical bubble collapse

The numerical implementation is verified by comparison with the analytical
solution of the Rayleigh-Plesset equation for spherical bubble dynamics. It
has influenced several works on the study of cavitation and bubble dynamics,
and has seen many adaptations over the last century. We present here only
salient features relevant to the current study and interested readers are
directed to comprehensive texts such as [7, 17].

We consider the case of an iso-thermal collapse of a spherical vaporous
bubble. A bubble of radius $R_{0}$ is initialized in an infinite domain of
liquid with a constant far-field pressure $p_{\infty}$. It is assumed that the
bubble contents are homogeneous and consist only of saturated vapor and no
non-condensable gases. It is also assumed that the pure liquid is
incompressible. In the absence of any non-condensable gases, the uniform
pressure $p_{B}(t)$ inside the bubble equals $p_{v}$. Fig. (2) shows the
schematic of the domain. Spherical symmetry is assumed and a one-dimensional
analysis is performed along the radial direction.

$\Omega^{\mathrm{f}}$$\Omega^{\mathrm{s}}$$\Gamma^{\mathrm{fs}}$collapsing
cavitiesresonating tip vortex cavityre-entrant jetshed cavitiesleading-edge
partial cavities$\&$ vortices$U_{\infty}$[A] Navier-Stokes Solver[C]
Cavitation SolverNon-
linear[D]$\boldsymbol{u}^{\mathrm{n+1}}_{\mathrm{(k+1)}},p^{\mathrm{n+1}}_{\mathrm{(k+1)}}$1\.
Interpolate solution2\. Solve3\. Correct solution4\. Update solution1\.
Interpolate solution2\. Solve3\. Correct solution4\. Update
solution$\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\alpha(\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k)}}-\boldsymbol{u}^{\mathrm{f},\mathrm{n}})$$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow
p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k)}}$$\Delta\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}$
and $\Delta p^{\mathrm{f},\mathrm{n+1}}$ in Eq.
(46)$\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}+\Delta\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}$$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow
p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}+\Delta
p^{\mathrm{f},\mathrm{n+1}}$$\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\dfrac{1}{\alpha}(\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}-\boldsymbol{u}^{\mathrm{f},\mathrm{n}})$$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow
p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}+\Delta
p^{\mathrm{f},\mathrm{n+1}}$$\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\phi^{\mathrm{f},\mathrm{n}}+\alpha(\phi^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k)}}-\phi^{\mathrm{f},\mathrm{n}})$$\Delta\phi^{\mathrm{f},\mathrm{n+\alpha}}$
in Eq.
(47)$\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}+\Delta\phi^{\mathrm{f},\mathrm{n+\alpha}}$$\phi^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow\phi^{\mathrm{f},\mathrm{n}}+\dfrac{1}{\alpha}(\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}-\phi^{\mathrm{f},\mathrm{n}})$iterations$\phi^{\mathrm{n+1}}_{\mathrm{(k+1)}}$[B]$r$$R(t)$$p(r,t)$$u(r,t)$$p_{B}(t)=p_{v}$$Liquid$$Liquid-
vapor\>interface$$Vapor$$r_{\infty},p_{\infty}$ Figure 2: Representative
schematic for iso-thermal collapse of a vaporous spherical bubble

Using mass conservation in the incompressible liquid outside the bubble, the
radial outward velocity $u(r,t)$ at a radial distance $r$ from the center of
the bubble can be shown to follow an inverse square law of the form

$u(r,t)\propto 1/r^{2}$ (53)

Assuming that the liquid at the interface is in thermal equilibrium with the
vapor at the saturation temperature, no mass transfer in the form of
evaporation or condensation occurs at the interface. Thus, at the interface
the velocity $u(R,t)=dR/dt$. Thus Eq. (53) can be written as

$u(r,t)=\frac{R^{2}}{r^{2}}\frac{dR}{dt}$ (54)

where $R(t)$ is the bubble radius. Momentum conservation in the $r$-direction
is given by

$-\frac{1}{\rho_{l}}\frac{\partial p}{\partial r}=\frac{\partial u}{\partial
t}+u\frac{\partial u}{\partial
r}-\nu_{l}\left[\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\frac{\partial u}{\partial r}\right)-\frac{2u}{r^{2}}\right]$
(55)

Substituting the expression for $u$ from Eq. (54) the following relation is
obtained

$-\frac{1}{\rho_{l}}\frac{\partial p}{\partial
r}=2\left(\frac{dR}{dt}\right)^{2}\left[\frac{R}{r^{2}}-\frac{R^{4}}{r^{5}}\right]+\frac{R^{2}}{r^{2}}\frac{d^{2}R}{dt^{2}}$
(56)

Assuming surface tension and viscous effects to be negligible, and the absence
of any mass transfer at the interface, $p_{r=R}=p_{B}(t)=p_{v}$. Using this
assumption, spatial integration of Eq. (56) between $r=R$ and
$r\rightarrow\infty$, gives the well-known Rayleigh equation [44]

$\frac{p_{v}-p_{\infty}}{\rho_{l}}=R\frac{d^{2}R}{dt^{2}}+\frac{3}{2}\left(\frac{dR}{dt}\right)^{2}$
(57)

Eq. (57) can be integrated to obtain the interface velocity during the
collapse

$\frac{dR}{dt}=-\sqrt{\frac{2}{3}\frac{p_{v}-p_{\infty}}{\rho_{l}}\left(1-\frac{R_{0}^{3}}{R^{3}}\right)}$
(58)

The interface acceleration is obtained as

$\frac{d^{2}R}{dt^{2}}=\frac{p_{v}-p_{\infty}}{\rho_{l}}\frac{R_{0}^{3}}{R^{4}}$
(59)

Using the kinematic relations in Eq. (58) and Eq. (59), Eq. (56) can be
integrated between $r$ (in the liquid outside the bubble) and
$r\rightarrow\infty$ to obtain the pressure $p(r,t)$ at the radial distance
$r$ from the centre of the bubble

$\frac{p(r,t)-p_{\infty}}{p_{\infty}-p_{v}}=-\frac{4}{3}\left(1-\frac{R_{0}^{3}}{R^{3}}\right)\left(\frac{R}{r}-\frac{1}{4}\frac{R^{4}}{r^{4}}\right)-\frac{R^{3}_{0}}{R^{2}r}$
(60)

In the absence of thermal effects and non-condensable gas content, we
encounter the special case of a vaporous bubble collapsing to zero volume. The
total time $t_{tc}$ taken by such a bubble to collapse from an initial radius
$R_{0}$ is presented by Rayleigh’s relation [44]:

$t_{tc}=0.915R_{0}\left(\frac{\rho_{l}}{p_{\infty}-p_{v}}\right)^{1/2}$ (61)

###  3.2 Numerical case setup

As we discussed in Section 1, the numerical modeling of cavitation requires
careful consideration of the cavitation model based on the spatial-temporal
scales of the particular flow configuration. We use cavitation model A for the
numerical study of this case of micro-scale bubble collapse because of its
origins in the Rayleigh-Plesset equation. In addition, model A has been
previously used to study this configuration in [19], where an implementation
based on the finite volume method was used. For comparison of the efficacy of
the present variational finite element implementation, we set up our numerical
case using similar geometrical parameters and model coefficients as in [19].

For the numerical study we consider a 3D spherical domain of
$R_{\infty}=0.5$\mathrm{m}$$, consisting of quiescent liquid. A spherical
vaporous bubble of radius $R_{0}=4\times 10^{-4}$\mathrm{m}$$ is initialized
in the centre of the domain. The spherical domain is discretized using
$840264$ hexahedral elements, with 49 nodes in the radial direction resolving
the initial bubble. The initial phase indicator and pressure inside the bubble
are set to $0.01$ and $p_{v}=2320$\mathrm{P}\mathrm{a}$$ respectively,
corresponding to the vapor phase. Outside the bubble, the phase indicator is
set to $1$. The pressure in the liquid is initialized according to Eq. (60).
Fig. (3) shows the initial conditions for the phase indicator. A far-field
pressure $P_{\infty}=1\times 10^{5}~{}$\mathrm{P}\mathrm{a}$$ is weakly
enforced over the outer surface boundary boundary of the domain using a
traction boundary condition. The densities of the pure liquid and vapor phases
are taken as $\rho_{l}=1000~{}$\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}$$ and
$\rho_{v}=0.01389~{}$\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}$$, giving a density
ratio $\dfrac{\rho_{l}}{\rho_{v}}\approx 72000$. The dynamic viscosities of
the two phases are set to
$\mu_{l}=0.001~{}$\mathrm{k}\mathrm{g}\,\mathrm{m}^{-1}\,\mathrm{s}^{-1}$$ and
$\mu_{v}=10^{-5}~{}$\mathrm{k}\mathrm{g}\,\mathrm{m}^{-1}\,\mathrm{s}^{-1}$$.
The model parameters $n_{0}$ and $d_{nuc}$ are assumed to be $10^{8}$ and
$10^{-4}~{}$\mathrm{m}$$ respectively. The condensation and evaporation
coefficients are set to $C_{c}=C_{v}=100$.

Figure 3: Sliced octant of domain showing the initial phase indicator values.
Inlay showing a close-up of the bubble at the centre of the domain. Inside the
bubble, the pressure is initialized to the vapor pressure $p_{v}$, with
$\phi^{\mathrm{f}}=0.01$ corresponding to the vapor phase. Outside the bubble,
the initial pressure has a continuous distribution according to Eq. (60) and a
constant $\phi^{\mathrm{f}}=1$ corresponding to the liquid phase.

###  3.3 Results and discussion

The results of the numerical study are first compared to the exact solution of
Eq. (58). Fig. (4(a)) shows the evolution of the bubble radius $R(t)$ with
time. In the numerical study, we consider the bubble radius to be the
effective radius of the volume of vapor in the domain.

$R=\left(\dfrac{3}{4\pi}\int_{\Omega^{\mathrm{f}}}\left(1-\phi^{\mathrm{f}}\right)\mathrm{d}\Omega^{\mathrm{f}}\right)^{1/3}$
(62)

Since no evaporation occurs outside the bubble, and no convection of vapor
occurs though the boundary at $R_{\infty}$
$\left(\dfrac{\partial\phi^{\mathrm{f}}}{\partial r}=0\right)$, the only vapor
content in the domain is in the bubble. This is also confirmed later in Fig.
(6(a)). Thus, Eq. (62) is seen to represent the bubble radius well. The bubble
radius and the solution time are non-dimensionalized using the initial bubble
radius $R_{0}$ and the Rayleigh collapse time $t_{tc}$ respectively. We
consider predictions obtained using four different values of the numerical
time-step $\Delta t$. It can be observed from Fig. (4(a)) that solutions
obtained using $\Delta t\leq 1\times 10^{-7}~{}s$ are able to capture the
evolution of the bubble radius well. To proceed, we sample the error between
the numerical and exact bubble radius at four instances during the collapse
process. These sampling instances correspond to the physical times
$0.405t_{tc}$, $0.810t_{tc}$, $0.945t_{tc}$ and $0.972t_{tc}$. We compute the
Root Mean Square Percentage Error (RMSPE) using the sampled results as

$\mathrm{RMSPE}=\sqrt{\frac{1}{n}\cdot\sum_{i=1}^{n}\Delta
R_{\mathrm{rel},i}^{2}}\cdot 100\%$ (63)

with

$\Delta R_{\mathrm{rel},i}=\frac{R_{num,i}}{R_{exact,i}}-1$ (64)

where $R_{num,i}$ is the numerically obtained bubble radius and $R_{exact,i}$
is the exact solution and $i$ denotes the index of the sampling time. Figure
(4(b)) shows the calculated RMSPE for different values of $\Delta t$. We
observe that the RMSPE is below $2.5\%$ for $\Delta t\leq 5\times
10^{-8}~{}$\mathrm{s}$$, indicating good agreement with the exact solution.
Thus, for the rest of this study we present results obtained with $\Delta
t=5\times 10^{-8}~{}$\mathrm{s}$$.

######  Remark 1.

The solution is observed to be convergent and stable across a range of time-
step sizes, even at large values of the order of $\Delta t=1\times
10^{-6}$\mathrm{s}$$. No numerical oscillations are seen in the solution
fields in the domain. [19] reported the presence of spurious pressure pulses
in the domain for $\Delta t>5\times 10^{-8}$\mathrm{s}$$. We note the efficacy
of the present implementation in suppressing these spurious oscillations.

(a)

(b)

Figure 4: Numerical assessment of spherical bubble collapse problem: (a)
comparison of bubble radius against analytical solution of the Rayleigh
equation (b) root mean square percentage error in bubble radius calculation at
different values of the numerical time-step size $\Delta t$

(a)

(b)

(c)

Figure 5: Spherical bubble collapse problem: (a) comparison of non-dimensional
pressure $\Pi$ with the analytical solution at different time steps, (b)
sensitivity of non-dimensional pressure $\Pi$ on the coefficients $C_{c}$ and
$C_{v}$, and (c) assessment of spherical symmetry.

(a)

(b)

Figure 6: Spherical bubble collapse problem: Distribution of
$\phi^{\mathrm{f}}$ and $\left(\dfrac{\partial\phi^{\mathrm{f}}}{\partial
r}\right)$ in the vicinity of the bubble surface

We next look at the predicted pressure field, which is another quantity of
interest in cavitating flows. It is well known that the cavity collapse
process can result in high pressures in the domain, several times the
magnitude of the collapse driving pressure $p_{\infty}$. Thus, for any
effective study of the fluid-structure interaction and noise effects in
cavitating flows, it is important that the cavitation solver is able to
accurately capture the pressure field. Figure (5(a)) shows the non-dimensional
pressure $\Pi=\dfrac{p^{\mathrm{f}}-p_{\infty}}{p_{\infty}-p_{v}}$ in the
domain at different times during the bubble collapse, compared to the
analytical solution of Eq. (60). The pressure $p^{\mathrm{f}}$ is taken in the
radial direction along the $z$-coordinate axis, with $r=0$ at the center of
the bubble. The numerically obtained pressures are observed to be in good
agreement with the analytical solution. The peak pressures are close to the
expected values, and an improvement in accuracy is obtained compared to
results presented in [19]. No spurious spikes or numerical oscillations are
observed in the pressure.

However, it is observed that at the final stages of the collapse process, the
pressure inside the bubble is higher than $p_{v}$. This deviates from the
assumption that the bubble pressure $p_{B}(t)$ equals $p_{v}$ at all times.
Increasing the coefficients $C_{c}$ and $C_{v}$ appears to resolve this. Fig.
(5(b)) compares the pressure profile at $t=0.945\tau_{R}$ for three values of
the evaporation and condensation coefficients. At $C_{c}=C_{v}=175$, the
pressure profile nearly matches the analytical solution. The pressure inside
the bubble is also observed to approach closer to the theoretical value of the
vapor pressure. On further increasing the $C_{c}=C_{v}$ to $250$, the vapor
pressure inside the bubble is recovered. However, the peak pressure is
observed to be under-predicted. Thus, a systematic tuning of the semi-
empirical coefficients is required, which is one of the limitations of such
phenomenological models. However, further tuning of these parameters is beyond
the scope of the current work.

Unlike in [19] where spherical symmetry was assumed in the numerical
simulation, in the current study we consider the full 3D spherical domain.
Numerical discretization and the weakly enforced far-field pressure boundary
condition can introduce asymmetry in the solution field. To determine the
ability of the solver to preserve the symmetrical nature of the solution, the
non-dimensional pressure $\Pi$ is plotted along the three cartesian
directions. Fig. (5(c)) shows the comparison of $\Pi$ (at $t=0.972t_{tc}$,
$C_{c}=C_{v}=100$) in the $x$, $y$ and $z$ directions as a function of the
distance from the centre of the bubble. It is observed that the 3D solver is
able to naturally preserve the symmetry of the collapsing bubble.

######  Remark 2.

We note the ability of the present implementation to accurately predict the
pressure in the domain, and the absence of spurious pressure oscillations
across the bubble interface. The vapor pressure inside the cavity is also
recovered, although it requires careful calibration of the model coefficients.
It is possible that an adaptive calibration (e.g., deep learning models [42,
9])) of these coefficients can aid in generalizing the model to multiple flow
configurations. This can be explored in future studies.

Next, we investigate the predicted values of the phase indicator
$\phi^{\mathrm{f}}$ in the domain. As discussed previously, presence of large
spatial gradients of $\phi^{\mathrm{f}}$ across the bubble interface can
result in unbounded numerical oscillations in the density field, and in turn
the pressure field. Fig. (6(a)) shows the distribution of $\phi^{\mathrm{f}}$
along the radial direction, while Fig. (6(b)) shows the gradient
$\left(\dfrac{\partial\phi^{\mathrm{f}}}{\partial r}\right)$ of
$\phi^{\mathrm{f}}$ with respect to the radial distance $r$. We observe no
oscillations in the solution in the vicinity of the bubble surface. The
solution is seen to be bounded within the range
$\phi^{\mathrm{f}}\in\left[0,1\right]$. The gradient
$\dfrac{\partial\phi^{\mathrm{f}}}{\partial r}$ is $\geq 0$ across the
interface, indicating monotone solutions. The location of the peak values of
the gradient are observed to be coincident with the exact value of the bubble
radius. Once again, a value of $C_{c}=C_{v}=175$ is observed to best represent
the expected distribution of $\phi^{\mathrm{f}}$ in the domain.

######  Remark 3.

We note the ability of the present implementation to recover monotone and
bounded solutions across the bubble interface. The absence of spurious
oscillations in $\phi^{\mathrm{f}}$ in the vicinity of the interface allows
the solver to be stable at large values of the density ratio
$\dfrac{\rho_{l}}{\rho_{v}}$. It is worth emphasizing, in the current study,
the density ratio is taken $\approx 72000$.

##  4 Turbulent cavitating flow over a hydrofoil: A validation study

In this section, we validate the proposed numerical method on the case of
turbulent cavitating flow over a hydrofoil. This is an often-encountered
scenario in marine propellers where fluid acceleration over the hydrofoil
surface can result in very low pressures and cavity inception near the blade
leading edge. For this study, we use cavitation model B because of its origins
in flows involving large bubble clusters. It is computationally less
expensive, and has been previously applied to the study of macro-scale
cavitation over hydrofoils. The turbulence model is validated first on non-
cavitating flow before proceeding to the case of cavitating flow. Fig. (7)
shows the general schematic of the computational domain used in the sections
to follow. $C$ is the hydrofoil chord length, $\alpha$ is the angle of attack
of the incoming flow, $H$ is the channel height and $\nu_{T}$ is the kinematic
turbulence viscosity. Specific details are given in the respective case
descriptions.

$\Omega^{\mathrm{f}}$$\Omega^{\mathrm{s}}$$\Gamma^{\mathrm{fs}}$collapsing
cavitiesresonating tip vortex cavityre-entrant jetshed cavitiesleading-edge
partial cavities$\&$ vortices$U_{\infty}$[A] Navier-Stokes Solver[C]
Cavitation SolverNon-
linear[D]$\boldsymbol{u}^{\mathrm{n+1}}_{\mathrm{(k+1)}},p^{\mathrm{n+1}}_{\mathrm{(k+1)}}$1\.
Interpolate solution2\. Solve3\. Correct solution4\. Update solution1\.
Interpolate solution2\. Solve3\. Correct solution4\. Update
solution$\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\alpha(\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k)}}-\boldsymbol{u}^{\mathrm{f},\mathrm{n}})$$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow
p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k)}}$$\Delta\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}$
and $\Delta p^{\mathrm{f},\mathrm{n+1}}$ in Eq.
(46)$\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}+\Delta\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}$$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow
p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}+\Delta
p^{\mathrm{f},\mathrm{n+1}}$$\boldsymbol{u}^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow\boldsymbol{u}^{\mathrm{f},\mathrm{n}}+\dfrac{1}{\alpha}(\boldsymbol{u}^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}-\boldsymbol{u}^{\mathrm{f},\mathrm{n}})$$p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow
p^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}+\Delta
p^{\mathrm{f},\mathrm{n+1}}$$\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\phi^{\mathrm{f},\mathrm{n}}+\alpha(\phi^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k)}}-\phi^{\mathrm{f},\mathrm{n}})$$\Delta\phi^{\mathrm{f},\mathrm{n+\alpha}}$
in Eq.
(47)$\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}\leftarrow\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}+\Delta\phi^{\mathrm{f},\mathrm{n+\alpha}}$$\phi^{\mathrm{f},\mathrm{n+1}}_{\mathrm{(k+1)}}\leftarrow\phi^{\mathrm{f},\mathrm{n}}+\dfrac{1}{\alpha}(\phi^{\mathrm{f},\mathrm{n+\alpha}}_{\mathrm{(k+1)}}-\phi^{\mathrm{f},\mathrm{n}})$iterations$\phi^{\mathrm{n+1}}_{\mathrm{(k+1)}}$[B]$r$$R(t)$$p(r,t)$$u(r,t)$$p_{B}(t)=p_{v}$$Liquid$$Liquid-
vapor\>interface$$Vapor$$r_{\infty},p_{\infty}$$\boldsymbol{u}^{\mathrm{f}}\cdot\boldsymbol{n}^{\mathrm{f}}=0,\boldsymbol{\sigma}\cdot\boldsymbol{n}^{\mathrm{f}}=0$$\boldsymbol{u}^{\mathrm{f}}=\left[U_{\infty},0,0\right]$$\boldsymbol{u}^{\mathrm{f}}\cdot\boldsymbol{n}^{\mathrm{f}}=0,\boldsymbol{\sigma}\cdot\boldsymbol{n}^{\mathrm{f}}=0$$Outflow$$\nabla\phi^{\mathrm{f}}\cdot\boldsymbol{n}^{\mathrm{f}}=0$$Inflow$$Slip$$Slip$$0.5C$$\alpha$$3.5C$$\phi^{\mathrm{f}}=1$$\boldsymbol{\sigma}\cdot\boldsymbol{n}^{\mathrm{f}}=0$$5.5C$$H$$y$$x$$\nabla\nu_{T}\cdot\boldsymbol{n}^{\mathrm{f}}=0$$No\textrm{-}slip$$\boldsymbol{u}^{\mathrm{f}}=0,\nu_{T}=0$
Figure 7: Representative computational domain and associated boundary
conditions for cavitating flow over hydrofoil

###  4.1 Turbulence model validation

In the current study, we employ a hybrid URANS-LES model to model turbulence.
We validate the turbulence model on non-cavitating flow over a hydrofoil
section. A NACA66 hydrofoil with chord length $C=0.15m$ and a span equal to
$0.3C$ is used in the study. The height of the channel $H=1.28C$. Fig. (8)
shows the computational grid. The grid consists of 76364 hexahedral elements
(eight-node bricks) resolving the cross-section. 30 nodes resolve the spanwise
direction. A target $y^{+}=yu_{\tau}/\nu=1$ was maintained in the
discretization of the hydrofoil boundary layer, where $y$ is the height of the
first node from the wall, $u_{\tau}$ is the friction velocity and $\nu$ is the
kinematic viscosity of the single phase liquid. The flow Reynolds number is
$8\times 10^{5}$, with a free-stream velocity
$U_{\infty}=5.333~{}$\mathrm{m}\,\mathrm{s}^{-1}$$. Liquid water at
$25^{\circ}C$ is taken as the working fluid, with a single-phase density of
$\rho_{l}=999.19~{}$\mathrm{k}\mathrm{g}\,\mathrm{m}^{-3}$$. A Dirichlet
velocity condition equal to $U_{\infty}$ is set at the inlet with a natural
traction-free outflow condition at the flow exit. A symmetric boundary
condition is used on the top and bottom surfaces with periodic conditions on
the spanwise surfaces. The time-step $\Delta t$ of the numerical simulation is
set to $1.407\times 10^{-4}$($T_{ref}/200$, where $T_{ref}=C/U_{\infty}$).

Figure 8: Computational mesh for NACA66 hydrofoil. Inlay showing mesh in the
vicinity of the hydrofoil

For our validation, the angle of attack $\alpha$ of the hydrofoil is varied
between $0^{\circ}-4^{\circ}$ for the non-cavitating flow condition, and the
time-averaged lift $(C_{L})$ and drag $(C_{D})$ coefficients were monitored.
Fig. (9) shows the comparison of $C_{L}$ and $C_{D}$ obtained from the
numerical simulation with the experimental results in [38]. The predicted
numerical results are observed to agree well with the experimental values in
the non-cavitating regime, and are within the uncertainties for $C_{L}(\Delta
C_{L}=0.012)$ and $C_{D}(\Delta C_{D}=0.002)$ reported in [38].

(a)

(b)

Figure 9: NACA66 hydrofoil problem: Predicted $C_{L}$ and $C_{D}$ in the non-
cavitating range compared with the experimental values from [38]

###  4.2 Cavitating flow over a hydrofoil

Figure 10: Computational mesh for NACA0012 hydrofoil. Inlay showing mesh in
the vicinity of the hydrofoil

To proceed further, the case of turbulent cavitating flow over a hydrofoil is
studied using our numerical solver. A NACA0012 hydrofoil with
$\alpha=1^{\circ}$ is consider with a Reynolds number of $2\times 10^{6}$. The
chord length $C=1$ with the height of the channel $H=6C$. Figure (10) shows
the computational grid used in the study. The domain is discretized with 90340
hexahedral elements and a 2D periodic boundary condition is applied in the
spanwise direction. A target $y^{+}$ equal to 1 is enforced at the hydrofoil
surface. The fluid domain is initialized with a liquid phase fraction of
$\phi^{\mathrm{f}}=1$. A freestream velocity of
$U_{\infty}=1$\mathrm{m}\,\mathrm{s}^{-1}$$ is applied at the inlet as a
dirichlet boundary condition. A traction-free outflow boundary condition is
used, weakly setting $p_{\infty}=0$. The phase fraction is set to $1$ at the
inlet, along with a Neumann boundary condition at the outflow. A density ratio
($\rho_{l}/\rho_{v}$) of 1000 is used in the study. The cavitation number of
the flow is defined as
$\sigma=\dfrac{p_{\infty}-p_{v}}{0.5\rho_{l}U_{\infty}^{2}}$, where
$p_{\infty}$ is the free-steam hydrostatic pressure. In the current study, the
cavitation number of the flow is set to $0.42$. The vapor pressure is set to
meet the cavitation number of the flow. A time-step of $\Delta t=1\times
10^{-4}$\mathrm{s}$$ ($t_{ref}/1000$) is utilized for the study.

Figure 11: Comparison of predicted pressure coefficient $C_{p}$ with results
presented in [49]

Figure (11) shows the result of the numerical study. The pressure coefficient
($C_{p}=\dfrac{p^{\mathrm{f}}-p_{\infty}}{0.5\rho_{l}U_{\infty}^{2}}$) on the
suction surface of the hydrofoil is compared against the numerical results of
[49] based on Kunz et al. [35] model. A good agreement can be seen between
these two numerical studies. The cavity pressure has been captured
consistently, which is important for the prediction of hydrodynamic loads on
the hydrofoil surface. A sharp gradient in the pressure over the hydrofoil
suction surface can be observed, corresponding to the cavity closure location.
Inlay shows the partial sheet cavity on the suction surface of the hydrofoil,
marked by an iso-contour of $\phi^{\mathrm{f}}=0.95$. Notably, it is observed
to stabilize over time to form a thin attached partial cavity on the hydrofoil
surface. A slight shift in the cavity closure location is observed which is
attributed to the difference in the cavitation and turbulence models used in
the two studies. In addition, the model coefficients used in the current study
are taken from[25] and [51] where the geometries studied were different. No
cavity separation and shedding are observed, and a fully attached turbulent
flow exists over the hydrofoil surface. The absence of a re-entrant jet is
consistent with the observations of [22] for thin cavities.

##  5 Application to fluid-structure interaction of caviating hydrofoil

In this closing section, we explore the ability of the proposed implementation
to model freely moving cavitating hydrofoil and to predict some key features
of cavitating flow over hydrofoils. We also test the compatibility of the
solver on configurations with moving solid boundaries. Before proceeding to
our fully-coupled cavitation and FSI demonstration, the primary objective of
this section is to evaluate the feasibility of the implementation for
turbulent cavitating flows over a stationary hydrofoil section. We consider
the same NACA0012 hydrofoil geometry as in Sec. (4.2) with the following
modifications. The height of the channel $H$ is reduced to $1.28C$. The
computational grid is also modified to a hybrid grid system comprising
hexahedral (8-node brick) and prism (6-node wedge) elements. This is to
prevent the formation of highly skewed elements (in fully hexahedral grids)
resulting from mesh deformations at high angles of attack. Figure (12) shows
the representative computational grid in the vicinity of the hydrofoil,
demonstrating deformation to $\alpha=10^{\circ}$. Cavitation model B is used
for the following studies. The values of the model coefficients and fluid
properties defined in Section (4.2) are used. The cavitation number $\sigma$
is increased to $1.2$ and the flow Reynolds number $Re$ is set to $1\times
10^{6}$, with $U_{\infty}=1$\mathrm{m}\,\mathrm{s}^{-1}$$. A traction-free
outflow is used, setting the pressure weakly to $0$ at the outflow boundary.

Figure 12: Hybrid computational grid for NACA0012 hydrofoil deformed to
$\alpha=10^{\circ}$

###  5.1 Stationary hydrofoil

Leading-edge cavitation over hydrofoils can behave in different ways based on
flow conditions defined by the flow Reynolds number $Re$, the cavitation
number $\sigma$ and the angle of attack $\alpha$ of the incoming flow. Under
certain combinations of these flow parameters, periodic cavity growth and
shedding can be observed. Integral to this cavity shedding process is the
formation of a periodic re-entrant jet along the hydrofoil surface. This jet
periodically flows along the hydrofoil surface from the cavity closure
location towards the leading edge, detaching the attached cavity. The
capturing of this flow phenomena is of interest to the study of propeller
vibration and noise because of two reasons. First, the periodic shedding of
the cavity leads to periodic fluctuations in the hydrodynamic loading on the
propeller blade. If the frequency of this loading is close to the natural
frequency of the blade, it can result in structural excitation - leading to
vibration and tonal noise emission. Second, the shed cavities exist in the
form of clouds of vaporous bubbles. These bubbles are prone to collapsing near
the trailing edge of the blade with localized high-amplitude water hammer
impacts. This can contribute to broadband underwater noise emission. We
investigate here the turbulent cavitating flow over a hydrofoil section at
$\alpha=10^{\circ}$.

(a) $t_{0}$

(b) $t_{0}+0.11t_{\infty}$

(c) $t_{0}+0.17t_{\infty}$

(d) $t_{0}+0.24t_{\infty}$

(e) $t_{0}+0.39t_{\infty}$

(f) $t_{0}+0.44t_{\infty}$

(g) $t_{0}+0.49t_{\infty}$

(h) $t_{0}+0.53t_{\infty}$

Figure 13: Contours of Z-vorticity (positive is clockwise, negative is anti-
clockwise) during one shedding cycle. Cavity (in red) marked by iso-contour of
$\phi^{\mathrm{f}}=0.95$. First figure marked by time $t_{0}$ for reference.
Subsequent figures marked in terms of $t_{0}$ and $t_{\infty}=C/U_{\infty}$

Figure 14: Vortex interaction and formation of re-entrant jet at trailing edge
of hydrofoil. Cavity (in red) marked at iso-contour of
$\phi^{\mathrm{f}}=0.95$.

Figure (13) shows the results of the study at different instances during one
cavity shedding cycle. An iso-contour of $\boldsymbol{\phi}^{\mathrm{f}}=0.95$
(in red) is used to represent the cavity surface. Also shown are the contours
of the vorticity in the direction $z$ out of the plane of the figure. Fig.
(13(a)) shows the inception of a leading-edge cavity. The cavity is observed
to grow, primarily collocated with a leading-edge vortex (LEV), to the extent
of the hydrofoil chord. In Fig. (13(d)) the clockwise rotating LEV interacts
with a counter-clockwise trailing edge vortex (TEV). We observe that the
interaction, along with a reverse pressure gradient, leads to the formation of
a re-entrant jet along the hydrofoil suction surface. This detaches the cavity
which is shed in the form of pockets (clouds) and is convected with the mean
flow. Fig. (14) shows the streamlines in the domain at the beginning of the
shedding process. The formation of the re-entrant jet can be observed to
originate at the intersection of the LEV and TEV.

######  Remark 4.

We note the ability of the present implementation to capture some select
physics of interest in turbulent cavitating flow over hydrofoils. Only
preliminary results are presented for demonstration. The accuracy of the
shedding frequency needs to be investigated, and can require careful
calibration of coefficients in the cavitation model and the modified turbulent
viscosity. Detailed investigation is beyond the scope of the current work, and
will be explored in future studies.

###  5.2 Pitching hydrofoil

Figure 15: Representative prescribed pitching motion

(a) $(+)5^{\circ}$

(h) $(-)5^{\circ}$

(b) $(+)10^{\circ}$

(g) $(-)10^{\circ}$

(c) $(+)12.5^{\circ}$

(f) $(-)12.5^{\circ}$

(d) $(+)14^{\circ}$

(e) $(-)14^{\circ}$

Figure 16: Contours of Z-vorticity (positive is clockwise, negative is anti-
clockwise) during one pitching cycle. Cavity (in red) marked by iso-contour of
$\phi^{\mathrm{f}}=0.95$. (+) marks the upward pitching stroke and (-) marks
the downward stroke.

We next investigate turbulent cavitating flow over the NACA0012 section
subject to a prescribed periodic pitching motion. By using our ALE-based FSI
framework, we linearly ramp the angle of attack of the hydrofoil between
$0^{\circ}$ and $15^{\circ}$. Fig. (15) shows the first four cycles of the
prescribed motion. The frequency $f_{pitch}$ of the motion is
$2~{}$\mathrm{H}\mathrm{z}$$. The rest of the study parameters are kept the
same as in Section (5.1). Figure (16) shows the cavities in the domain marked
by the iso-contour $\phi^{\mathrm{f}}=0.95$. Also plotted are the contours of
z-vorticity. At the high pitching frequency
$f_{pitch}=2~{}$\mathrm{H}\mathrm{z}$$ and the associated hydrofoil
acceleration, cavities are observed to originate on both surfaces because of
the low pressures during the pitching motion. The primary cavity generation is
at the hydrofoil leading edge. The cavity shedding frequency is low compared
to the pitching frequency and the hydrofoil suction surface is seen to be
perennially covered by cavities that are continuously being shed and convected
with the mean flow. The collocation between the vortices and the cavities can
be discerned. Detailed investigations on the influence of the pitching
frequency on the cavity and vortex shedding frequency are reserved for future
work.

We note the compatibility of the cavitation and the flow solvers with
structural deformation. The versatility provided by the partitioned coupling
was exploited to couple the additional solvers for turbulence modeling and ALE
mesh update. The numerical solution obtained is stable and $2-3$ non-linear
predictor-corrector iterations are seen to give a converged solution at each
time-level. This sets the stage for large-scale fully coupled FSI studies in
cavitating flows, and will be of interest for future work.

##  6 Conclusion

A robust and accurate variational finite element formulation for the numerical
study of cavitating flows has been presented for stationary and moving
hydrofoils. We introduced novel stabilized linearizations of two cavitation
transport equation models based on the two-phase homogeneous mixture theory.
The numerical implementation has been employed to study two cavitating flow
configurations with vastly different temporal and spatial scales. An initial
verification study on the micro-scale collapse of a spherical vaporous bubble
has been shown to maintain stability across a range of time steps. Accurate
solutions of the pressure field were obtained, devoid of spurious numerical
oscillations. The solution of the phase indicator is demonstrated to be
bounded, and the solver is numerically stable at large density ratios. The
implementation has been validated on fully turbulent cavitating flow over a
hydrofoil with good agreement with previous numerical and experimental
studies. We also explored the ability of the implementation to predict select
characteristic features of macro-scale cavitating flows, including re-entrant
jets, periodic cavity shedding and cavity-vortex interaction. Further, we
examined the versatility of the implementation to be coupled with solvers for
studying fluid-structure interaction, demonstrating the case of a pitching
hydrofoil. In future work, the authors plan to extend the implementation to
study fully-coupled FSI studies in cavitating flows. One potential application
is cavitating flow-induced vibrations taking into account the hydroelastic
response of structures. Another potential application is the study of material
erosion resulting from cavitation bubble collapse using fully-Eulerian fluid-
solid formulations.

## Acknowledgements

The authors would like to acknowledge the Natural Sciences and Engineering
Research Council of Canada (NSERC) for the funding. This research was enabled
in part through computational resources and services provided by (WestGrid)
(https://www.westgrid.ca/), Compute Canada (www.computecanada.ca) and the
Advanced Research Computing facility at the University of British Columbia.

## References

  * [1] MPI: A message-passing interface version 3.1.  (www.mpi-formum.org).  Technical report, 2015.
  * Akcabay and Young [2014] Deniz Tolga Akcabay and Yin Lu Young.  Influence of cavitation on the hydroelastic stability of hydrofoils.  _Journal of Fluids and Structures_ , 49:170–185, 2014.
  * Arakeri et al. [1988] VH Arakeri, H Higuchi, and REA Arndt.  A model for predicting tip vortex cavitation characteristics.  _Journal of Fluids Engineering_ , 110(2):190–193, 1988.
  * Arndt et al. [2015] R Arndt, P Pennings, J Bosschers, and T Van Terwisga.  The singing vortex.  _Interface focus_ , 5(5):20150025, 2015.
  * Bailey et al. [2003] Michael R Bailey, Robin O Cleveland, Tim Colonius, Lawrence A Crum, Andrew P Evan, James E Lingeman, James A McAteer, Oleg A Sapozhnikov, and JC Williams.  Cavitation in shock wave lithotripsy: the critical role of bubble activity in stone breakage and kidney trauma.  In _IEEE Symposium on Ultrasonics, 2003_ , volume 1, pages 724–727. IEEE, 2003.
  * Bayram and Korobenko [2020] A Bayram and A Korobenko.  Variational multiscale framework for cavitating flows.  _Computational Mechanics_ , 66(1):49–67, 2020\.
  * Brennen [2013] Christopher Earls Brennen.  _Cavitation and Bubble Dynamics_.  Cambridge University Press, 2013.  doi: 10.1017/CBO9781107338760.
  * Brooks and Hughes [1982] Alexander N Brooks and Thomas JR Hughes.  Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations.  _Computer methods in applied mechanics and engineering_ , 32(1-3):199–259, 1982.
  * Bukka et al. [2021] Sandeep Reddy Bukka, Rachit Gupta, Allan Ross Magee, and Rajeev Kumar Jaiman.  Assessment of unsteady flow predictions using hybrid deep learning based reduced-order models.  _Physics of Fluids_ , 33(1):013601, 2021.
  * Carlton [2018] John Carlton.  _Marine propellers and propulsion_.  Butterworth-Heinemann, 2018.
  * Chahine et al. [2016] G. L. Chahine, A Kapahi, J.-K. Choi, and C.-T. Hsiao.  Modeling of surface cleaning by cavitation bubble dynamics and collapse.  _Ultrasonics Sonochemistry_ , 29:528–549, 2016.
  * Chahine et al. [1983] Georges L Chahine, Andrew F Conn, Virgil E Johnson Jr, et al.  Cleaning and cutting with self-resonating pulsed water jets.  In _2nd US water jet conference_. Citeseer, 1983.
  * Chen and Heister [1996] Yongliang Chen and SD Heister.  Modeling hydrodynamic nonequilibrium in cavitating flows.  _Journal of Fluids Engineering_ , 118(1):172–178, 1996.
  * Chung and Hulbert [1993] Jintai Chung and GM Hulbert.  A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method.  _Journal of Applied Mechanics_ , 1993.
  * Coutier-Delgosha et al. [2003] O Coutier-Delgosha, JL Reboud, and Y Delannoy.  Numerical simulation of the unsteady behaviour of cavitating flows.  _International journal for numerical methods in fluids_ , 42(5):527–548, 2003.
  * Deshpande et al. [1994] Manish Deshpande, Jinzhang Feng, and Charles L Merkle.  Cavity flow predictions based on the euler equations.  _Journal of Fluids Engineering_ , 1994.
  * Franc and Michel [2006] Jean-Pierre Franc and Jean-Marie Michel.  _Fundamentals of cavitation_ , volume 76.  Springer science & Business media, 2006.
  * Franca and Frey [1992] Leopoldo P Franca and Sérgio L Frey.  Stabilized finite element methods: Ii. the incompressible navier-stokes equations.  _Computer Methods in Applied Mechanics and Engineering_ , 99(2-3):209–233, 1992.
  * Ghahramani et al. [2019] Ebrahim Ghahramani, Mohammad Hossein Arabnejad, and Rickard E Bensow.  A comparative study between numerical methods in simulation of cavitating bubbles.  _International Journal of Multiphase Flow_ , 111:339–359, 2019.
  * Gnanaskandan and Mahesh [2015] A Gnanaskandan and Krishnan Mahesh.  A numerical method to simulate turbulent cavitating flows.  _International Journal of Multiphase Flow_ , 70:22–34, 2015\.
  * Goncalves et al. [2010] Eric Goncalves, Maxime Champagnac, and Regiane Fortes Patella.  Comparison of numerical solvers for cavitating flows.  _International Journal of Computational Fluid Dynamics_ , 24(6):201–216, 2010.
  * Gopalan and Katz [2000] Shridhar Gopalan and Joseph Katz.  Flow structure and modeling issues in the closure region of attached cavitation.  _Physics of fluids_ , 12(4):895–911, 2000.
  * Harari and Hughes [1992] Isaac Harari and Thomas JR Hughes.  What are c and h?: Inequalities for the analysis and design of finite element methods.  _Computer Methods in Applied Mechanics and Engineering_ , 97(2):157–192, 1992.
  * Hsu et al. [2010] M-C Hsu, Yuri Bazilevs, Victor M Calo, Tayfun E Tezduyar, and Thomas JR Hughes.  Improving stability of stabilized and multiscale formulations in flow simulations at small time steps.  _Computer Methods in Applied Mechanics and Engineering_ , 199(13-16):828–840, 2010.
  * Huang et al. [2013] Biao Huang, Antoine Ducoin, and Yin Lu Young.  Physical and numerical investigation of cavitating flows around a pitching hydrofoil.  _Physics of Fluids_ , 25(10):102109, 2013.
  * Hughes and Wells [2005] Thomas JR Hughes and Garth N Wells.  Conservation properties for the galerkin and stabilised forms of the advection–diffusion and incompressible navier–stokes equations.  _Computer methods in applied mechanics and engineering_ , 194(9-11):1141–1159, 2005.
  * Husseini et al. [2005] Ghaleb A Husseini, Mario A Diaz De La Rosa, Eric S Richardson, Douglas A Christensen, and William G Pitt.  The role of cavitation in acoustically activated drug delivery.  _Journal of Controlled Release_ , 107(2):253–261, 2005.
  * Jansen et al. [2000] Kenneth E Jansen, Christian H Whiting, and Gregory M Hulbert.  A generalized-$\alpha$ method for integrating the filtered navier–stokes equations with a stabilized finite element method.  _Computer methods in applied mechanics and engineering_ , 190(3-4):305–319, 2000.
  * Ji et al. [2015] B Ji, XW Luo, Roger EA Arndt, Xiaoxing Peng, and Yulin Wu.  Large eddy simulation and theoretical investigations of the transient cavitating vortical flow structure around a naca66 hydrofoil.  _International Journal of Multiphase Flow_ , 68:121–134, 2015.
  * Jiang et al. [2018] Jingwei Jiang, Haopeng Cai, Cheng Ma, Zhengfang Qian, Ke Chen, and Peng Wu.  A ship propeller design methodology of multi-objective optimization considering fluid–structure interaction.  _Engineering Applications of Computational Fluid Mechanics_ , 12(1):28–40, 2018.
  * Joshi and Jaiman [2017a] Vaibhav Joshi and Rajeev K Jaiman.  A positivity preserving variational method for multi-dimensional convection–diffusion–reaction equation.  _Journal of Computational Physics_ , 339:247–284, 2017a.
  * Joshi and Jaiman [2017b] Vaibhav Joshi and Rajeev K Jaiman.  A variationally bounded scheme for delayed detached eddy simulation: Application to vortex-induced vibration of offshore riser.  _Computers & fluids_, 157:84–111, 2017b.
  * Joshi and Jaiman [2018] Vaibhav Joshi and Rajeev K Jaiman.  A positivity preserving and conservative variational scheme for phase-field modeling of two-phase flows.  _Journal of Computational Physics_ , 360:137–166, 2018\.
  * Kerr et al. [1940] W Kerr, JF Shannon, and RN Arnold.  The problems of the singing propeller.  _Proceedings of the Institution of Mechanical Engineers_ , 144(1):54–90, 1940.
  * Kunz et al. [1999] Robert F Kunz, David A Boger, Thomas S Chyczewski, D Stinebring, H Gibeling, and T Govindan.  Multi-phase cfd analysis of natural and ventilated cavitation about submerged bodies.  In _Proceedings of the 3rd ASME-JSME Joint Fluids Engineering Conference_ , 1999.
  * Lee et al. [2017] Abe H Lee, Robert L Campbell, Brent A Craven, and Stephen A Hambric.  Fluid–structure interaction simulation of vortex-induced vibration of a flexible hydrofoil.  _Journal of Vibration and Acoustics_ , 139(4), 2017.
  * Lee et al. [2014] Hyoungsuk Lee, Min-Churl Song, Jung-Chun Suh, and Bong-Jun Chang.  Hydro-elastic analysis of marine propellers based on a bem-fem coupled fsi algorithm.  _International journal of naval architecture and ocean engineering_ , 6(3):562–577, 2014.
  * Leroux et al. [2004] Jean-Baptiste Leroux, Jacques André Astolfi, and Jean Yves Billard.  An experimental study of unsteady partial cavitation.  _Journal of Fluids Engineering_ , 126(1):94–101, 2004.
  * Maines and Arndt [1997] B Maines and Roger EA Arndt.  The case of the singing vortex.  _Journal of Fluids Engineering_ , 1997.
  * Maljaars et al. [2018] Pieter Maljaars, Mirek Kaminski, and Henk Den Besten.  Boundary element modelling aspects for the hydro-elastic analysis of flexible marine propellers.  _Journal of Marine Science and Engineering_ , 6(2):67, 2018.
  * Merkle [1998] Charles L Merkle.  Computational modelling of the dynamics of sheet cavitation.  In _Proc. of the 3rd Int. Symp. on Cavitation, Grenoble, France, 1998_ , 1998.
  * Miyanawala and Jaiman [2017] Tharindu P Miyanawala and Rajeev K Jaiman.  An efficient deep learning technique for the navier-stokes equations: Application to unsteady wake flow dynamics.  _arXiv preprint arXiv:1710.09099_ , 2017.
  * Pishchalnikov et al. [2003] Yuriy A Pishchalnikov, Oleg A Sapozhnikov, Michael R Bailey, James C Williams Jr, Robin O Cleveland, Tim Colonius, Lawrence A Crum, Andrew P Evan, and James A McAteer.  Cavitation bubble cluster activity in the breakage of kidney stones by lithotripter shockwaves.  _Journal of endourology_ , 17(7):435–446, 2003\.
  * Rayleigh [1917] Lord Rayleigh.  Viii. on the pressure developed in a liquid during the collapse of a spherical cavity.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 34(200):94–98, 1917.
  * Ross and Kuperman [1989] Donald Ross and WA Kuperman.  Mechanics of underwater noise, 1989.  URL https://doi.org/10.1121/1.398685.
  * Saad and Schultz [1986] Youcef Saad and Martin H Schultz.  Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  _SIAM Journal on scientific and statistical computing_ , 7(3):856–869, 1986.
  * Schnerr and Sauer [2001] Günter H Schnerr and Jürgen Sauer.  Physical and numerical modeling of unsteady cavitation dynamics.  In _Fourth international conference on multiphase flow_ , volume 1. ICMF New Orleans, 2001.
  * Schnerr et al. [2008] Günter H Schnerr, Ismail H Sezal, and Steffen J Schmidt.  Numerical investigation of three-dimensional cloud cavitation with special emphasis on collapse induced shock dynamics.  _Physics of Fluids_ , 20(4):040703, 2008.
  * Senocak and Shyy [2001] Inanc Senocak and Wei Shyy.  Numerical simulation of turbulent flows with sheet cavitation.  2001\.  URL http://resolver.caltech.edu/cav2001:sessionA7.002.
  * Senocak and Shyy [2002] Inanc Senocak and Wei Shyy.  A pressure-based method for turbulent cavitating flow computations.  _Journal of Computational Physics_ , 176(2):363–383, 2002.
  * Senocak and Shyy [2004] Inanc Senocak and Wei Shyy.  Interfacial dynamics-based modelling of turbulent cavitating flows, part-1: Model development and steady-state computations.  _International journal for numerical methods in fluids_ , 44(9):975–995, 2004.
  * Sezal [2009] Ismail Hakki Sezal.  _Compressible dynamics of cavitating 3-D multi-phase flows_.  PhD thesis, Technische Universität München, 2009.
  * Shakib et al. [1991] Farzin Shakib, Thomas JR Hughes, and Zdeněk Johan.  A new finite element formulation for computational fluid dynamics: X. the compressible euler and navier-stokes equations.  _Computer Methods in Applied Mechanics and Engineering_ , 89(1-3):141–219, 1991.
  * Singhal et al. [2002] 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, 2002.
  * Song et al. [2004] WD Song, Hong MH, Lukyanchuk B, and Chong TC.  Laser-induced cavitation bubbles for cleaning of solid surfaces.  _Journal of Applied Physics_ , 95:2952––2956, 2004.
  * Stride and Coussios [2019] Eleanor Stride and Constantin Coussios.  Nucleation, mapping and control of cavitation for drug delivery.  _Nature Reviews Physics_ , 1(8):495–509, 2019\.
  * Tezduyar et al. [1992] Tayfun E Tezduyar, Sanjay Mittal, SE Ray, and R Shih.  Incompressible flow computations with stabilized bilinear and linear equal-order-interpolation velocity-pressure elements.  _Computer Methods in Applied Mechanics and Engineering_ , 95(2):221–242, 1992.
  * Van Oossanen [1974] Pieter Van Oossanen.  _Calculation of performance and cavitation characteristics of propellers including effects on non-uniform flow and viscosity_.  PhD thesis, Delft University of Technology, 1974.  URL http://resolver.tudelft.nl/uuid:daef4d65-e0cc-4796-88e0-6c6b2d2a54ac.
  * Wu et al. [2018] Qin Wu, Biao Huang, Guoyu Wang, and Shuliang Cao.  The transient characteristics of cloud cavitating flow over a flexible hydrofoil.  _International Journal of Multiphase Flow_ , 99:162–173, 2018.
  * Zijlstra [2011] AG Zijlstra.  _Acoustic surface cavitation_.  PhD thesis, University of Twente, 2011.
  * Zwart et al. [2004] Philip J Zwart, Andrew G Gerber, Thabet Belamri, et al.  A two-phase flow model for predicting cavitation dynamics.  In _Fifth international conference on multiphase flow, Yokohama, Japan_ , volume 152, 2004.