1. Introduction
Noise pollution due to human activity (i.e. anthropogenic noise) poses serious threats to the marine ecosystem (Duarte et al. Reference Duarte2021). One of the main sources of anthrophony is shipping, which generates noise through the operation of propellers, hull vibrations and onboard machinery (Smith & Rigby Reference Smith and Rigby2022). A significant proportion of the propeller noise is due to the cavitation phenomenon occurring in various forms, such as sheet and tip vortex cavitation (TVC). In most cases, tip vortex cavitation is the first form of cavitation appearing on ship propellers (Zhang et al. Reference Zhang, Zhang, Peng, Wang and Shao2015), which arises due to the low pressure within the vortex trailing from the tip of the propeller blades.
The tip vortex forms on lifting surfaces of finite length due to the pressure difference between the pressure side and the suction side, which disappears at the tip. The pressure gradient on the pressure and suction sides leads to opposite spanwise velocity components on these surfaces, which generates a vortex that trails from the tip (Platzer & Souders Reference Platzer and Souders1979). The wake sheet shed from the blade rolls up into the tip vortex as the flow moves downstream, further strengthening the tip vortex (Batchelor Reference Batchelor1964). The pressure within the vortex core may drop below the vapour pressure, leading to a cavitating tip vortex. Tip vortex cavitation encompasses a complex dynamical interaction among vortical motion, cavitation and turbulence (Arndt Reference Arndt2002). This complex phenomenon can be observed in a relatively simple hydrofoil section that produces a strong swirling flow with pressure reduction and the cavitation process. This paper is motivated by the need to provide an improved understanding of the interplay of cavitation and vortex dynamics.
Understanding the physics of TVC – including the tip vortex roll-up and formation, TVC inception, tip vortex cavity dynamics, and the contribution of TVC to the underwater radiated noise (URN) level in various operating conditions – is crucial for the development of effective TVC mitigation strategies. To study the tip vortex flow field, one needs to consider numerous aspects of multiscale spatial-temporal dynamics, even in non-cavitating (wetted) conditions. During the formation of the tip vortex, the tip vortex/boundary layer interactions strongly affect the tip vortex flow field (Maines & Arndt Reference Maines and Arndt1997). Downstream of the tip, the influence of the roll-up of the wake in the tip vortex and the complex system of streamwise vortices shed from the tip become dynamically relevant (Devenport et al. Reference Devenport, Rife, Liapis and Follin1996). Further downstream, viscous decay and the introduction of vortex instabilities play an important role in the breakdown of the tip vortex (Ganesh, Schot & Ceccio Reference Ganesh, Schot and Ceccio2014).
The occurrence of cavitation further complicates the phenomenon of the tip vortex. These complexities coupled with the high three-dimensionality and anisotropy of the flow with extreme gradients (Asnaghi, Svennberg & Bensow Reference Asnaghi, Svennberg and Bensow2020), the impact of turbulent fluctuations on the cavitation inception (Arndt & Keller Reference Arndt and Keller1992), and tip vortex wandering, which gives rise to the need for special treatments when analysing the acquired data (Devenport et al. Reference Devenport, Rife, Liapis and Follin1996), make the investigation of tip vortex cavitating flows a complex endeavour. Despite all of these difficulties associated with the investigation of cavitating tip vortex flow, various researchers have studied the TVC phenomenon using theoretical, experimental and numerical methods focusing on the contributions of TVC to URN.
Recently, there has been a growing interest in studying the contributions of TVC to URN. Specifically, a significant focus has been placed on investigating the oscillations of the interface between the liquid phase and the vapour filament forming within the tip vortex, referred to as the cavity surface in the rest of this paper. Different modes of oscillation, namely, the breathing mode, serpentine (displacement) mode and double-helical mode, are observed on the tip vortex cavity interface. These modes of cavity oscillations are depicted in figure 1. The cavity oscillations have been found to affect the radiated noise (Bosschers Reference Bosschers2018b; Wang et al. Reference Wang, Bai, Cheng, Ji and Peng2023). An analytical solution was derived by Thomson (Reference Thomson1880) for the dispersion relation of excitations on a hollow columnar vortex representing a vortex cavity. For the purpose of acoustic analysis, Morozov (Reference Morozov1974) extended this solution to compressible flows. Bosschers (Reference Bosschers2018b) further improved the solution by including viscous effects and constant axial velocity. Experimental evidence for this analytical solution was further explored by Pennings et al. (Reference Pennings, Bosschers, Westerweel and van Terwisga2015a) for the tip vortex cavitating flow over a stationary NACA66(2)-415 hydrofoil with an elliptical planform. Furthermore, different modes of TVC were observed to induce a range of instabilities in the tip vortex, thus altering its dynamics (Ye, Wang & Shao Reference Ye, Wang and Shao2023). Of particular interest in our study is to investigate the development of cavity surface oscillatory behaviour as the tip vortex cavity travels downstream, to gain insight into the oscillation dynamics of different regions of the cavity and their contribution to pressure fluctuations in the surroundings.
Figure 1. Three main tip vortex cavity oscillation modes: (a) breathing, (b) displacement and (c) double-helical.
Due to the limitations of experimental methods in cavitating flow investigations, numerical studies can reveal more details about the flow and cavity characteristics in TVC analysis. However, obtaining sufficiently accurate results has been a recent development because of the challenges associated with simulating tip vortex cavitating flows. Large eddy simulation (LES) was shown to outperform Reynolds-averaged Navier–Stokes models for simulation of the tip vortex flow field by Asnaghi, Bensow & Svennberg (Reference Asnaghi, Bensow and Svennberg2017), and Asnaghi et al. (Reference Asnaghi, Svennberg and Bensow2020) successfully replicated the flow field of the tip vortex in wetted conditions observed in the experiments of Pennings, Westerweel & Terwisga (Reference Pennings, Westerweel and Terwisga2015b) using implicit LES. In cavitating conditions, modelled based on homogeneous mixture theory, their results regarding the azimuthal velocity distribution and the diameter of the cavitating tip vortex also agreed with the experimental data. They also proposed mesh resolution requirements for accurate simulation of tip vortex cavitating flows in terms of the number of nodes across the vortex core. To account for the effect of non-condensable gas bubbles on TVC, Cheng et al. (Reference Cheng, Long, Ji, Peng and Farhat2021) developed an Euler–Lagrangian cavitation model in which the gas bubbles’ motion was simulated through a Lagrangian approach, and their impact on cavitation was taken into account by modifying the cavitation model employed in the Eulerian-based simulation of the cavitating flow based on the local partial pressure of the non-condensable gas.
To the best of the authors’ knowledge, the first numerical study focusing on tip vortex cavity surface oscillation was carried out by Klapwijk et al. (Reference Klapwijk, Lloyd, Vaz, van den Boogaard and van Terwisga2022) using partially averaged Navier–Stokes and delayed detached eddy simulation. Their results included some of the cavity surface oscillations, but the breathing mode oscillations were not captured. Furthermore, it is evident from their results of the cavity morphology that the tip vortex cavity was over-dissipative and could not capture a major proportion of the surface oscillations. The double-helical mode of cavity surface oscillations was captured in the simulations of Wang et al. (Reference Wang, Bai, Qian, Cheng and Ji2022). They investigated the effect of roll-up of secondary vortices into the tip vortex flow, which was shown to be responsible for the enhancement of tip vortex instability. Recently, the cavity surface oscillations of the breathing mode were successfully captured in the numerical study by Wang et al. (Reference Wang, Bai, Cheng, Ji and Peng2023). The dominance of the hydroacoustic contribution of the breathing mode was also demonstrated in their work, with the tonal frequency and the centre frequency of the broadband hump observed in the noise spectrum being attributed to the breathing mode at zero streamwise wavenumber and the tip vortex cavity resonance. However, their study includes the results for only one cavitation number. Moreover, they considered only the developed region of the tip vortex cavity and not the formation region, i.e. where the wake roll-up into the tip vortex is not complete. In addition, their method for extracting the breathing mode oscillations from the cavity surface oscillation data requires improvement due to the centre of the cavity not being well defined in tip vortex cavitating flows.
In this study, using LES and the homogeneous mixture model, which has been shown to be able to predict the surface oscillations of the tip vortex cavity in previous studies (Klapwijk et al. Reference Klapwijk, Lloyd, Vaz, van den Boogaard and van Terwisga2022; Wang et al. Reference Wang, Bai, Qian, Cheng and Ji2022, Reference Wang, Bai, Cheng, Ji and Peng2023), the tip vortex cavitating flow over a NACA66(2)-415 stationary hydrofoil is investigated numerically via our in-house three-dimensional (3-D) variational finite element solver. In order to address the lack of a general estimation method for mesh resolution requirements in TVC simulation, a length scale is defined based on the radial pressure gradient within the vortex core, which takes into account the strength of the tip vortex as well as its core radius. This length scale is then employed to non-dimensionalize the mesh resolution in the tip vortex region to propose a generalizable method for the estimation of the required mesh resolution in the simulation of tip vortex flow in various configurations, rather than merely considering the number of grid nodes across the vortex core radius. After the assessment of the numerical results with the experimental and analytical data, a comprehensive investigation of the breathing mode of cavity surface oscillations is carried out for representative cavitation numbers. To extract the breathing mode oscillations, we define an effective radius based on the cross-sectional area of the tip vortex cavity that captures the cavity volume variations. Different regions as the tip vortex cavity progresses away from the tip are examined using the temporally averaged effective radius. The spatially averaged effective radius is analysed to discern the temporal variations of the vapour volume of the entire tip vortex cavity. The spatial-temporal characteristics of local breathing mode oscillations are further investigated to understand the characteristics of this oscillatory mode and its development along the cavity in different cavitation numbers. Proper orthogonal decomposition (POD) is employed to further examine the streamwise variations of the characteristics of the breathing mode. The pressure fluctuations are probed at different locations within the domain to investigate how the contributions of TVC to the pressure fluctuations vary as a function of cavitation numbers. To understand the correlations between the breathing oscillations and the pressure fluctuations, the spectra of the pressure fluctuations are compared to the breathing mode oscillation spectra within different regions of the tip vortex cavity. These findings may serve as a foundation for the development of mitigation strategies aimed at addressing tip vortex cavitation noise.
The remainder of the paper is structured as follows. In § 2, the computational framework employed in this work is described. The problem setting and the grid generation strategy are presented in § 3. The results obtained from the simulations are covered in § 4. While the focus in § 4.1 is placed on the non-cavitating conditions, § 4.2 provides an overview of the oscillation dynamics of the tip vortex cavity. A systematic investigation of the breathing mode of oscillation is presented in § 4.3. The pressure fluctuations in different cavitation numbers and their correlation with breathing mode oscillations are discussed in § 4.4. The key findings and concluding remarks are summarized in § 5.
2. Computational framework
In this work, our in-house high-fidelity finite element flow solver was employed for simulating the tip vortex cavitating flow. The mathematical modelling and numerical formulation are briefly presented here. For further details, the reader is referred to the original papers on the implementation of the models used herein (Jaiman, Guan & Miyanawala Reference Jaiman, Guan and Miyanawala2016; Kashyap & Jaiman Reference Kashyap and Jaiman2021).
2.1. Multiphase flow modelling
In this work, multiphase flow is treated as a continuous homogeneous mixture of liquid and vapour phases. A phase indicator 
$\phi ^f(\boldsymbol {x},t)$ determines the phase fraction of the liquid phase at any coordinate and time in the physical domain of the working fluid 
$\varOmega ^f(\boldsymbol {x},t)$, where 
$\boldsymbol {x}$ and 
$t$ are the spatial and temporal coordinates within the domain 
$\varOmega ^f(\boldsymbol {x},t)$ with the boundary 
$\varGamma^{{\kern-0.7pt}f}$. The mixture density (
$\rho ^f$) and the mixture dynamic viscosity (
$\mu ^f$) are considered to be linear combinations of the density and viscosity of the liquid and vapour phases, and are calculated as
$$\begin{gather} \rho^f = \rho_l \phi^f + \rho_v (1-\phi^f) , \end{gather}$$
$$\begin{gather}\mu^f = \mu_l \phi^f + \mu_v (1-\phi^f) , \end{gather}$$
where the properties with 
$l$ and 
$v$ subscripts represent the properties of the pure liquid and pure vapour phases, respectively. In order to obtain the value of the phase indicator 
$\phi ^f$ within the domain, a scalar transport equation is solved, which is given by
\begin{equation} \frac{\partial \phi^f}{\partial t} + \phi^f \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi^f = \frac{\dot{m}}{\rho_l} \quad \textrm{on}\ (\boldsymbol{x},t) \in \varOmega^f , \end{equation}
where 
$\boldsymbol {u}$ is the fluid velocity at each coordinate 
$(\boldsymbol {x}, t)$, and 
$\dot {m}$ is the finite mass transfer rate from the vapour phase to the liquid phase occurring due to cavitation.
2.2. Cavitation modelling
In (2.3), the source term 
$\dot {m}$ representing the mass transfer rate from the vapour phase to the liquid phase, which includes the effects of both condensation and evaporation, is calculated using the model proposed by Schnerr & Sauer (Reference Schnerr and Sauer2001):
\begin{align} \dot{m} & = \frac{3 \rho_l \rho_v}{\rho^f R_B} \sqrt{\frac{2}{3 \rho_l\,|p^f - p_v|}} \left[ C_c \phi^f (1-\phi^f) \max(p^f - p_v, 0) \right.\nonumber\\ &\quad \left.{}+ C_v \phi^f (1+\phi_{nuc} - \phi^f) \min(p^f - p_v, 0) \right] , \end{align}
where 
$p^f (\boldsymbol {x}, t)$ is the working fluid pressure, and 
$p_v$ is the saturation pressure of the working fluid. Also, 
$C_c$ and 
$C_v$ are the condensation and evaporation coefficients, which were added to the model in later numerical implementations (Cazzoli et al. Reference Cazzoli, Falfari, Bianchi, Forte and Catellani2016; Ghahramani, Arabnejad & Bensow Reference Ghahramani, Arabnejad and Bensow2019) in order to modify the condensation and evaporation behaviour with respect to specific configurations.
In this model, the initial nuclei within the flow are assumed to be of equal radius and homogeneously distributed within the domain. In (2.4), 
$\phi _{nuc}$ is the initial phase fraction of the bubble nuclei, which is related to the initial volume of the nuclei per unit volume as
\begin{equation} \phi_{nuc} = \frac{v_{nuc}}{1 + v_{nuc}} , \end{equation}
where 
$v_{nuc}$ is the initial volume of the nuclei per unit volume, which is related to the nuclei diameter (
$d_{nuc}$) and the number of nuclei per unit volume (
$n_0$) as 
$v_{nuc} = {{\rm \pi} n_0 d_{nuc}^3}/{6}$. In (2.4), 
$R_B$ is the equivalent radius of the vapour volume assumed to be in the form of one bubble, and is related to other parameters of the flow using the equation
\begin{equation} R_B = \left( \frac{3}{4 {\rm \pi}n_0}\,\frac{1 + \phi_{nuc} - \phi^f}{\phi ^ f} \right) ^ {1/3} . \end{equation}In the framework employed in this work, the effect of diffusion of dissolved air from the liquid phase into the tip vortex cavity is assumed to be negligible. It was found by Amini et al. (Reference Amini, Reclari, Sano, Dreyer and Farhat2019) that the diffusion of air into the tip vortex cavity, which causes hysteresis in the inception and desinence of TVC, is dependent mainly on the formation of a laminar separation bubble at the tip. Their results demonstrated that such a laminar separation bubble acts as a shelter for the tip vortex cavity and promotes the diffusion of air into the tip vortex cavity. They eliminated hysteresis caused by air diffusion into the tip vortex cavity through artificial roughening of the blade surface, which indicates that in the absence of the laminar separation bubble at the tip, the impact of diffusion of air into the tip vortex cavity is negligible. Since, as explained in § 3, the hydrofoil employed in the current work is a NACA 6 series with low adverse pressure gradients over its surface (Asnaghi et al. Reference Asnaghi, Svennberg and Bensow2020), no laminar separation bubble occurs at the tip of this hydrofoil, thus making it plausible to assume that the diffusion of dissolved air into the tip vortex cavity is negligible.
2.3. Fluid conservation of mass and momentum
The spatially filtered Navier–Stokes equations for an incompressible flow are
$$\begin{gather} \rho^f\,\frac{\partial \bar{\boldsymbol{u}}} {\partial t} + \rho^f \bar{\boldsymbol{u}} \boldsymbol{\cdot}\boldsymbol{\nabla} \bar{\boldsymbol{u}} = \boldsymbol{\nabla}\boldsymbol{\cdot} \bar{\boldsymbol{\sigma}}^f + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}^{sgs} +\boldsymbol{b}^f \quad \textrm{on} \ (\boldsymbol{x},t) \in \varOmega^f , \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \bar{\boldsymbol{u}} = 0 \quad \textrm{on} \ \varOmega^f , \end{gather}$$
where the overbar denotes the spatially filtered quantity, 
$\boldsymbol {b}^f$ denotes the body force exerted on the fluid, and the Cauchy and subgrid stresses are represented by 
$\bar {\boldsymbol {\sigma }}^f$ and 
$\boldsymbol {\sigma }^{sgs}$, respectively. The Cauchy stress for a Newtonian fluid can be written as
\begin{equation} \bar{\boldsymbol{\sigma}}^f ={-}\bar{p} \boldsymbol{\mathsf{I}} + \mu^f \left( \boldsymbol{\nabla} \bar{\boldsymbol{u}} + \left( \boldsymbol{\nabla} \bar{\boldsymbol{u}} \right) ^{\rm T} \right) , \end{equation}
where 
$\bar {p}$ represents the filtered fluid pressure. The subgrid stress 
$\boldsymbol{\sigma}^{sgs}$ is the additional stress term appearing in (2.7) due to the filtering applied in LES. In order to complete the set of equations, a model is required for the subgrid stress tensor.
2.4. Subgrid-scale modelling
Dynamic subgrid-scale (SGS) model is utilized in this work for modelling the SGS stress 
$\boldsymbol {\sigma }^{sgs}$. In this model, a filter, denoted by an overbar, is applied to the Navier–Stokes equations (2.7) and (2.8), which, as explained in § 2.3, gives rise to a new stress term, the SGS stress 
${\sigma }^{sgs} _ {ij} = \overline {u_i^f u_j^f} - \bar {u}_i^f \bar {u}_j^f$, appearing in the momentum equation. This stress term corresponds to the stresses originating from the scales of motion smaller than the spatial resolution 
$\varDelta$. As mentioned in § 2.3, since the SGS quantities are unknown, a model for the SGS stress is required for the closure of the problem.
The nonlinear SGS stress (Gatski & Speziale Reference Gatski and Speziale1993; Wang & Bergstrom Reference Wang and Bergstrom2005) is obtained using the equation
\begin{align} \sigma_{ij}^{sgs} -
\frac{\delta_{ij}}{3}\,\sigma_{kk}^{sgs} &\approx{-}2 \mu_t
\bar{{\mathsf{S}}}_{ij} - C_{NL} 6 \mu_t^2 /
\sigma_{kk}^{sgs}\nonumber\\ &\quad \times \left(
\bar{{\mathsf{S}}}_{ik} \bar{\varOmega}_{kj} + \bar{{\mathsf{S}}}_{jk}
\bar{\varOmega}_{ki} - 2 \bar{{\mathsf{S}}}_{ik} \bar{{\mathsf{S}}}_{kj} +
\frac{2}{3}\,\bar{{\mathsf{S}}}_{nk} \bar{{\mathsf{S}}}_{kn} \delta _ {ij}
\right) ,
\end{align}
where 
$\mu _t$ is the dynamic eddy viscosity, and
\begin{equation} \bar{{\mathsf{S}}}_{ij} \equiv \frac{1}{2} \left(\frac{\partial \bar{u}_i^f}{\partial x_j} + \frac{\partial \bar{u}_j^f}{\partial x_i}\right)\quad {\rm and}\quad \bar{\varOmega}_{ij} \equiv \frac{1}{2} \left(\frac{\partial \bar{u}_i^f}{\partial x_j} - \frac{\partial \bar{u}_j^f}{\partial x_i}\right) \end{equation}
are the resolved strain rate and rate of rotation, respectively. The eddy viscosity 
$\mu _t$ in this equation is related to the mesh resolution and the resolved strain tensor as
\begin{equation} \mu_t = \rho^f (C_s \bar{\varDelta})^2\,|\bar{\boldsymbol{\mathsf{S}}}| , \end{equation}
where 
$|\bar {\boldsymbol{\mathsf{S}}}|$ is the norm of the resolved strain rate tensor, calculated as 
$|\bar {\boldsymbol{\mathsf{S}}}| = (2 \bar {{\mathsf{S}}}_{ij} \bar {{\mathsf{S}}}_{ij})^{1/2}$. In order for LES to yield reasonable results close to the wall, a simple algebraic eddy viscosity model (Balaras, Benocci & Piomelli Reference Balaras, Benocci and Piomelli1996; Piomelli & Balaras Reference Piomelli and Balaras2002) is employed instead of (2.12) close to the wall as 
${\mu _t}/{\mu ^f} = \kappa y^+_w (1 - {\rm e} ^ {-y_w^+} / A)$, where 
$y^+_w = y_w u_\tau / \nu$ is the non-dimensional distance to the wall in wall units, with the friction velocity and the model coefficient being denoted by 
$u_\tau$ and 
$\kappa$, respectively, and 
$A$ having constant value 19.
In the dynamic SGS model (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991), a test filter is defined in addition to the grid filter. The grid filter, denoted by an overbar, has scale dimension 
$\varDelta$, and the test filter, denoted by the operation 
$\widehat {(\;)}$, is any coarser-level filter with scale dimension 
$\hat {\varDelta }$. The SGS stress at the test level 
${\mathsf{T}}_{ij} = \widehat {\overline {u_i u_j}} - \hat {\bar {u}}_i \hat {\bar {u}}_j$ is obtained by test-filtering the equations of motion. The stress tensor 
$\mathscr {L}_{ij}$, which is the stress associated with the length scales between the test filter size and the grid filter size,
\begin{equation} \mathscr{L}_{ij} = \widehat{\bar{u}_i \bar{u}_j} - \hat{\bar{u}}_i \hat{\bar{u}}_j, \end{equation}
can be related to 
${\mathsf{T}}_{ij}$ and 
$\sigma ^{sgs}_{ij}$ as
\begin{equation} \mathscr{L}_{ij} = {\mathsf{T}}_{ij} - \sigma^{sgs}_{ij} . \end{equation}
Using the Smagorinsky eddy-viscosity model for the unknown stresses 
$\sigma ^{sgs}_{ij}$ and 
${\mathsf{T}}_{ij}$ yields
\begin{equation} \mathscr{L}_{ij} ={-} 2
C_s^2 \left( \hat{\varDelta}^2\,|\hat{\bar{\boldsymbol{\mathsf{S}}}}|\,
\hat{\bar{{\mathsf{S}}}}_{ij} - \varDelta^2\,| \bar{\boldsymbol{\mathsf{S}}} |\,\bar{{\mathsf{S}}}_{ij}\right),
\end{equation}
and the Smagorinsky constant is related to a dynamic coefficient 
$C_{ds}$ as 
$C_s = \sqrt {C_{ds}}$ for 
$C_{ds}>0$, and 
$C_s=0$ otherwise. The dynamic coefficient 
$C_{ds}$ (Lilly Reference Lilly1992) is calculated as
\begin{equation} C_{ds} =
\frac{1}{2}\,\frac{\langle {\mathsf{M}}_{ij} \mathscr{L}_{ij}
\rangle}{\langle {\mathsf{M}}_{lk} {\mathsf{M}}_{lk} \rangle} ,
\end{equation}
where 
${\mathsf{M}}_{ij} = \hat {\varDelta } ^ 2\,|\hat {\bar {\boldsymbol{\mathsf{S}}}}|\,\hat {\bar {{\mathsf{S}}}}_{ij} - \bar {\varDelta } ^ 2\, \widehat {|\bar {\boldsymbol{\mathsf{S}}}| \bar {{\mathsf{S}}}_{ij}}$, and 
$\langle \cdot \rangle$ denotes some type of smoothing process such as averaging.
3. Problem set-up
In this work, the cavitating flow over a stationary NACA66(2)-415 hydrofoil with an elliptical planform is considered. The geometry of the problem and the working conditions are presented in § 3.1, and the strategy employed for grid generation is explained in § 3.2, where the derivation of a pressure-gradient-based length scale is presented and is employed for non-dimensionalization of the mesh size within the tip vortex flow region for mesh requirements specification.
3.1. Geometry and working conditions
The case investigated numerically in this study involves the flow over a stationary NACA66(2)-415 hydrofoil, which is placed within a tunnel with a square cross-section. The geometry and dimensions of the domain and the hydrofoil, depicted in figure 2, are chosen to be identical to those of the experiments by Pennings et al. (Reference Pennings, Westerweel and Terwisga2015b) since their experimental data are considered herein for the validation of our numerical framework. Due to the manufacturing constraints encountered in the experimental study (Pennings et al. Reference Pennings, Westerweel and Terwisga2015b), the foil employed in this investigation is truncated at the trailing edge. The root chord length of the foil after truncation is 
$c_0 = 0.1256 \ \mathrm {m}$, and the half-span of the foil is 
$h = 0.150 \mathrm {m}$. As shown in figure 2, the inflow and outflow boundaries are located 
$5c_0$ and 
$10c_0$ in front of and behind the hydrofoil, respectively. The tip of the foil is placed at the centre of the test section at a distance 
$0.150 \ \mathrm {m}$ from the lateral walls. The origin of the coordinate system is considered at the tip in the simulations, and 
$x$, 
$y$ and 
$z$ directions correspond to the streamwise, transverse and spanwise directions, respectively, as shown in figure 2.
Figure 2. Schematic of the computational domain and boundary conditions.
The inflow velocity and the Reynolds number are kept constant at 
$U_\infty =6.8 \ \mathrm {m} \mathrm {s}^{-1}$ and 
$Re= {U_\infty c_0}/{\nu } = 8.95 \times 10^5$, respectively, in all cases, similar to the experiments. The simulations are carried out at angle of attack 
$\alpha =9^\circ$. We consider both non-cavitating (wetted) and cavitating conditions for our physical investigation. The cavitation number (
$\sigma$) is defined based on the outflow pressure as
\begin{equation} \sigma = \frac{p_\infty - p_{sat}}{\dfrac{1}{2} \rho U_{\infty} ^ 2} , \end{equation}
where 
$p_{sat}$ and 
$p_\infty$ are the saturation and the outflow pressure, respectively. To explore the oscillatory dynamics of the cavity surface, we consider three cavitation numbers, 
$\sigma =1.2$, 
$1.7$ and 
$2.6$. The primary emphasis is on the case with cavitation number 
$\sigma =1.7$ due to the occurrence of moderate cavitation in this scenario. The inflow velocity is set to 
$\boldsymbol {u}_{in} = (U_\infty, 0, 0)$. The outflow boundary condition is weakly (Neumann) set as constant pressure. The bottom wall of the domain and the hydrofoil surface are specified as no-slip walls, and other domain walls are of symmetry (slip wall) boundary condition.
To achieve statistically stationary conditions in a shorter computation time, the initial condition was obtained from simulations with a coarser grid and a larger time step, and the results were projected onto the desired grid. Using this initial condition, in the first step, the simulations were carried out for 11 
$t_{ref}$, i.e. the flow characteristic time scale 
${c_0}/{U_\infty }$, with time step 
$\Delta t = 5.41 \times 10^{-3} \, t_{ref}$ in wetted conditions. Subsequently, the time step was decreased to 
$\Delta t = 5.41 \times 10^{-4} \, t_{ref}$, which was found to be small enough for accurate results. Cavitation was started after 16 
$t_{ref}$ of flow time, and cavitating flow simulations were carried out for 38 
$t_{ref}$. The results obtained for the last 22 
$t_{ref}$ of the flow time were used in this work for post-processing and analysis. Due to the computation time limitations, the cases 
$\sigma = 2.6$ and 
$1.2$ were simulated for 30 
$t_{ref}$ and 22 
$t_{ref}$, respectively, following the onset of cavitation. For both of these cases, the last 11 
$t_{ref}$ of the flow time were employed for the analysis.
3.2. Grid generation strategy
It is well known that due to the high gradients and small scales in tip vortex cavitating flows, adequate mesh resolution in simulating such flows is of critical significance due to the high susceptibility of such flows to numerical diffusion (Asnaghi et al. Reference Asnaghi, Svennberg and Bensow2020). Therefore, an adequately fine mesh is required within the tip vortex flow region. Due to the various resolution requirements in different regions of cavitating tip vortex flow, the grid needs to be meticulously designed to satisfy mesh requirements while avoiding excessive computational costs.
The cross-sectional views of the grid distribution employed in this work in the streamwise and transverse directions are illustrated in figure 3. Triangular mesh elements are generated on the hydrofoil surface and extruded normal to the hydrofoil surface, generating prism layers as the boundary layer mesh. Outside the boundary layer mesh, tetrahedral elements are utilized. A refinement region is created enclosing the hydrofoil and the region where flow dynamics is present. A refinement box is added in the wake region to capture the effects of wake dynamics on the tip vortex roll-up process. A cylindrical refinement region with radius 
$r = 30 \ \mathrm {mm}$ is created, starting from a streamwise position slightly upstream of the tip and extending 
$2 c_0$ downstream of the tip. A smaller refinement region along the path of the tip vortex is created enclosing the tip vortex core to capture the extreme gradients in the tip vortex flow. This region is created by sweeping the tip vortex trajectory with a circle of radius 
$r = 5 \ \mathrm {mm}$, which is obtained by preliminary simulations with a coarse cylindrical refinement.
Figure 3. Grid distribution in the streamwise and cross-flow directions: (a) streamwise distribution, (b) in-plane distribution.
The mesh resolution within the tip vortex region is one of the major factors determining the modelling accuracy of TVC. Many researchers have proposed various mesh requirements for simulating TVC (Ahmad, Proctor & Perry Reference Ahmad, Proctor and Perry2013; Asnaghi et al. Reference Asnaghi, Svennberg and Bensow2020). All of these recommendations are based on the number of grid points across the tip vortex diameter; however, merely considering the number of grid points across the tip vortex diameter does not take into account all the influencing parameters, and would not necessarily lead to sufficient resolution. Other characteristics of the tip vortex flow, such as the circulation (
$\varGamma$), should also be taken into account. There might exist tip vortex flows with the same vortex core diameters but different strengths, and the vortex flow with higher strength, indeed, requires finer mesh resolution. Furthermore, the vortex core radius might not be known prior to the simulations, which gives rise to the need for a method to estimate the mesh resolution requirements for tip vortex flow simulation based on a priori known parameters. Therefore, a dimensionless mesh resolution within the tip vortex flow region (
$\Delta r^*$) is defined herein to account for the effect of the tip vortex strength as well as its core radius.
In the simulation of tip vortex cavitating flows, capturing the radial gradient of pressure 
${\partial p}/{\partial r}$ is of utmost importance. Therefore, the maximum radial pressure gradient together with the fluid density 
$\rho$ and kinematic viscosity 
$\nu$ is used to define a length scale for the tip vortex flow as
\begin{equation} \mathcal{L}_\varGamma = \left( \frac{\rho \nu ^ 2}{\left(\dfrac{\partial p}{\partial r}\right)_{max}} \right) ^{{1}/{3}} . \end{equation}
Since the maximum gradient of pressure is not known before the simulations and experiments, it needs to be related to a priori known properties of the flow, such as the Reynolds number and the lift coefficient (
$C_L$), which is known from the experimental data or can be obtained easily using low-cost numerical analyses when experimental data are not available. To relate the maximum radial gradient of pressure to the known parameters of the flow, the tip vortex flow field is assumed to be axisymmetric, and the radial velocity is assumed to be zero. In this case, the radial momentum equation simplifies to 
${\partial p}/{\partial r} = \rho ({u_\theta ^2}/{r})$, where 
$r$ is the radial distance from the vortex centre, and 
$u_\theta$ is the azimuthal velocity of the flow. The Rankine vortex model with circulation 
$\varGamma$ and core radius 
$r_v$ is employed to obtain the azimuthal velocity. The Rankine vortex model assumes a solid-body-like rotation within an inner core (
$r< r_v$), thereby having a radial velocity distribution 
$u_\theta (r) = {\varGamma r}/{2 {\rm \pi}r_v^2}$. The vortex flow outside the rigid core is assumed to be an irrotational vortex, i.e. 
$u_\theta (r) = {\varGamma }/{2 {\rm \pi}r}$.
The maximum radial pressure gradient in a Rankine vortex occurs on the vortex core boundary, i.e. 
$r=r_v$, where the azimuthal velocity in the Rankine vortex model is 
$u_\theta = {\varGamma }/{2 {\rm \pi}r_v}$, the substitution of which in the simplified radial momentum equation yields
\begin{equation} \left( \frac{\partial p}{\partial r} \right) _{max} = \frac{\rho \varGamma^2}{4 {\rm \pi}^2 r_v^3} \end{equation}
for the maximum radial pressure gradient. Using the maximum radial pressure gradient obtained from (3.3) in (3.2), and rearranging based on the vortex Reynolds number (
$Re_\varGamma = \varGamma / \nu$), yields
\begin{equation} \mathcal{L}_\varGamma = (2{\rm \pi})^{2/3} r_v\,Re_\varGamma ^ {-{2}/{3}} . \end{equation}
Based on (3.4), the length scale 
$\mathcal {L}_\varGamma$ scales linearly with the vortex core radius 
$r_v$, and scales inversely with 
$Re_\varGamma ^ {2/3}$, indicating that as the vortex core radius decreases or the vortex strength (represented by the vortex Reynolds number 
$Re_\varGamma$) increases, the length scale 
$\mathcal {L}_\varGamma$ decreases, thus indicating the need for a finer mesh to be employed in the tip vortex region.
The circulation of the tip vortex 
$\varGamma$ and the tip vortex core radius 
$r_v$ are unknown before the numerical analysis or experiments. Therefore, these quantities also need to be estimated for the calculation of the length scale 
$\mathcal {L}_\varGamma$. Fruman & Dugue (Reference Fruman and Dugue1992) assumed that the local tip vortex circulation and vortex core radius are related to the wing mid-span circulation 
$\varGamma _0$ and the boundary layer thickness 
$\delta$, respectively. The mid-span circulation of the hydrofoil 
$\varGamma _0$ is calculated using the Kutta–Jukowski theorem as 
$\varGamma _ 0 = \frac {1}{2} C_L U_\infty c_0$. Astolfi, Fruman & Billard (Reference Astolfi, Fruman and Billard1999) stated that the radius of the vortex core is directly related to the thickness of the turbulent boundary layer (
$\delta$) of a flat plate with length one blade root chord length 
$c_0$, which can be related to the Reynolds number using the equation 
$\delta = 0.37 c_0\,{Re}^{-0.2}$. Pennings et al. (Reference Pennings, Bosschers, Westerweel and van Terwisga2015a) obtained the proportionality coefficients 
$\varGamma /\varGamma _0$ and 
$r_v/\delta$ to be in the ranges 
$0.44<\varGamma /\varGamma _0<0.49$ and 
$0.56< r_v/\delta <0.61$. In the calculations of Astolfi et al. (Reference Astolfi, Fruman and Billard1999), these values were in the ranges 
$0.5<\varGamma /\varGamma _0<0.6$ and 
$0.8< r_v/\delta <1.0$. Since the goal is to find an estimate for the required mesh size within the tip vortex core, the turbulent boundary layer thickness and the mid-span wing circulation can be substituted for the tip vortex core radius and the tip vortex circulation in (3.4), respectively. Therefore, the coefficients can be omitted as they merely affect the scaling of the results. The pressure-gradient-based length scale 
$\mathcal {L}_\varGamma$ is described in terms of the a priori known flow parameters as
\begin{equation} \mathcal{L}_\varGamma = \left( \frac{16 {\rm \pi}^2 }{C_L^2\,{Re}^{2.6}} \right) ^ {{1}/{3}} c_0 . \end{equation}
The length scale 
$\mathcal {L}_\varGamma$ is employed to non-dimensionalize the mesh element size within the tip vortex region 
$\Delta r$ as 
$\Delta r ^ * = {\Delta r}/{\mathcal {L}_\varGamma }$. In this work, simulations are carried out using different values for the non-dimensional mesh resolution 
$\Delta r ^ *$ within the tip vortex region, and the mesh requirements are proposed accordingly in terms of the non-dimensional mesh resolution. The mesh resolutions within the tip vortex flow region and their corresponding non-dimensional values are summarized in table 1, the results of which are presented in § 4.1.
Table 1. Simulation cases for the mesh resolution effect study.
4. Results and discussion
In this section, the numerical results obtained using our finite element solver are presented. In § 4.1, our findings from non-cavitating tip vortex flow simulations are compared with the experimental data (Pennings et al. Reference Pennings, Westerweel and Terwisga2015b), and the effect of mesh resolution on the accuracy of the simulations is investigated systematically. The tip vortex cavitating flow simulation results are examined in § 4.2, and evaluated with the semi-analytical solution by Bosschers (Reference Bosschers2018b) as well as the experimental measurements of Pennings et al. (Reference Pennings, Westerweel and Terwisga2015b). The breathing mode of oscillation is further explored for different cavitation numbers in § 4.3, and the pressure fluctuations within the domain are analysed in § 4.4 to determine the correlations between the breathing mode oscillations of the tip vortex cavity and the pressure fluctuations within the domain.
4.1. Non-cavitating flow
To study the effect of mesh resolution on the wetted tip vortex flow simulation, and to validate the capability of the presented numerical framework in simulating such flows, the axial (
$x$ direction) and in-plane (parallel with the 
$yz$-plane) velocity distribution at three different downstream locations (
$x/c_0 = 0.5, 0.75, 1.14$) are compared qualitatively with the experimental data, and quantitative comparison is carried out for the time-averaged lift coefficient and the in-plane velocity distribution at the downstream location 
$x/c_0 = 1.14$. Due to a direct relationship between the lift coefficient and the tip vortex flow field, both the lift coefficient and the velocity profiles require accurate predictions for the modelling of tip vortex cavitating flow. The simulations obtain a time-averaged lift coefficient 
${C_L = 0.68}$, which is in agreement with the value reported in Pennings et al. (Reference Pennings, Westerweel and Terwisga2015b), i.e. 
${C_L \approx 0.66 \pm 0.02}$.
The velocity data obtained from the simulations should be processed before comparing with the experimental data. First, the centre of the vortex is identified at every time step assuming that the minimum pressure within the vortex core occurs at the centre of the vortex. In the next step, the results obtained from 3-D simulations are projected onto a two-dimensional (2-D) cross-section of the flow orthogonal to the streamwise direction (parallel with the 
$yz$-plane). The mesh projected onto the 2-D plane is illustrated in figure 4(a). In the next step, a new polar coordinate system is defined with the vortex core centre as the origin, and the data are interpolated from the 2-D plane of interest onto the polar mesh, depicted in figure 4(b), around the tip vortex core centre using the cubic spline interpolation method. This process is carried out for different flow times, then the quantities obtained for various flow times on the vortex-based polar coordinate system are averaged after aligning the tip vortex centre obtained for all flow times used in the averaging. Utilizing this method of post-processing eliminates the effect of vortex wandering on time averaging. For the quantitative comparison of azimuthal velocity at the downstream location 
$x/c_0 = 1.14$ with the experimental measurements, the in-plane velocity values are spatially averaged on 
$90^\circ$ arcs between the two perpendicular black lines illustrated in 5(b), the same as in the averaging method used in Pennings et al. (Reference Pennings, Westerweel and Terwisga2015b) since the asymmetry of the tip vortex flow field at this location necessitates the utilization of the same averaging method for consistency.
Figure 4. Projected 2-D triangular mesh and polar mesh used during the data processing for comparison with the experimental data. (a) Triangular mesh on a 2-D plane orthogonal to the free-stream flow direction obtained from interpolation from the 3-D mesh. (b) Polar grid points generated around the tip vortex centre for interpolation of the data from the 2-D cross-section of the flow.
Figure 5. Comparison of non-dimensional velocity distribution within the tip vortex region at the downstream location 
$x/c_0 = 1.14$ obtained from simulations with different mesh sizes with the experimental measurements. (a) Non-dimensional axial velocity. (b) Non-dimensional in-plane velocity.
As a quantitative comparison, the time- and contour-averaged azimuthal velocity profiles predicted by the simulations employing various mesh sizes within the tip vortex region at the cross-section 
$x/c_0 = 1.14$ are shown in figure 6. Pennings et al. (Reference Pennings, Westerweel and Terwisga2015b) employed two averaging methods to obtain the azimuthal velocity profile, namely sum-of-correlation (SOC) and conditional averaging, the latter of which is more similar to the method of averaging utilized in the present study. Their results obtained using both methods are plotted in figure 6. This plot shows that accurate capturing of the velocity field inside the tip vortex requires extremely high mesh resolution. The use of inadequate resolution can lead to excessive dissipation and underprediction of the maximum azimuthal velocity and its gradient. Moreover, employing adequately small grid elements within the tip vortex region is necessary for the correct prediction of the tip vortex core radius. The results presented in figure 6 indicate that the results obtained from the simulation with 
$\Delta r ^ * = 46$ are excessively dissipated compared to the experimental data. Refinement of the grid to 
$\Delta r ^ * = 23$ led to the correct overall capture of the velocity distribution, but the tip vortex flow field, in this case, is slightly dissipated. The utilization of a grid with 
$\Delta r ^ * = 18.5$ provides an accurate prediction of the tip vortex flow field.
Figure 6. Comparison of the normalized azimuthal velocity profile at 
$x/c_0 = 1.14$ obtained from simulations using various grid sizes with experimental data.
Based on the qualitative and quantitative comparisons of the predicted axial and azimuthal velocity fields presented herein, it can be concluded that the non-dimensional grid size 
$\Delta r ^ *$ is required to be below 20 for a reliable capture of tip vortex flow dynamics. This guideline can serve as a more generalizable method for estimating the required mesh resolution based on a priori known parameters for future works on tip vortex flow simulation compared to previously proposed mesh resolution requirements. Furthermore, the agreement between the results obtained from the simulation using a non-dimensional mesh resolution 
$\Delta r ^ * = 18.5$ in the tip vortex region reveals the capability of the developed finite-element-based numerical framework in simulating non-cavitating tip vortex flows.
4.2. Simulation of TVC oscillatory dynamics
Similar to the experimental observations (Pennings et al. Reference Pennings, Bosschers, Westerweel and van Terwisga2015a), the tip vortex cavity observed in the simulations carried out in this work exhibits complex surface oscillations, as shown in figure 7. As explained in § 1, these surface oscillations have been found to be one of the sources of TVC noise. Therefore, correctly capturing these oscillations in the simulations, and then investigating the physics of these oscillations and deformations, especially the breathing mode of oscillation, is crucial for understanding the contributions of the TVC to the URN. In this subsection, the effect of mesh resolution on cavity surface oscillations observed in the simulations is discussed in § 4.2.1, which is followed by an overview of the oscillations of the tip vortex cavity surface for the case 
$\sigma =1.7$ obtained from the simulation with 
$\Delta r^* = 18.5$.
Figure 7. Tip vortex cavity surface observed in simulations (viewed from the pressure side of the hydrofoil) at cavitation number 
$\sigma = 1.7$ with (a) 
$\Delta r^* = 46$, (b) 
$\Delta r^* = 23$, and (c) 
$\Delta r^* = 18.5$.
4.2.1. Effect of mesh resolution on TVC oscillatory dynamics
The impact of mesh resolution on the cavity surface oscillations is investigated via instantaneous tip vortex cavity shapes for various mesh resolutions. The tip vortex cavities depicted in figure 7 are extracted from the simulation results using the iso-surface of 
$\alpha = 0.9$, where 
$\alpha$ is the volume fraction of liquid water. As is evident in this figure, the tip vortex cavity obtained from simulations using the coarsest mesh resolution in this study, i.e. 
$\Delta r^* = 46$, does not exhibit physically consistent surface oscillations and lacks the double-helical mode specifically. Although the breathing and displacement modes are present in this case, they occur at excessively small amplitudes.
As the mesh resolution is refined to 
$\Delta r^* = 23$, surface deformations of the cavity captured in the simulations are intensified, manifesting a more realistic behaviour than the case with 
$\Delta r^* = 46$. Figure 7(b) demonstrates that the tip vortex cavity obtained using mesh resolution 
$\Delta r^* = 23$ possesses the stationary double-helical shape in the roll-up region of the tip vortex similar to the stationary double-helical mode of deformation in the time-averaged tip vortex cavity structure obtained from simulations using 
$\Delta r^* = 18.5$ displayed in figure 8(a). Furthermore, various modes of travelling surface oscillations are observed on the tip vortex cavity obtained from the simulation with 
$\Delta r^* = 23$ along the cavity.
Figure 8. Time-averaged cavitating tip vortex flow field at 
$\sigma = 1.7$. (a) Time-averaged tip vortex cavity obtained from simulations with 
$\Delta r^* = 18.5$. (b) Iso-surfaces of 
$(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla }) u ( {U_\infty ^2}/{c_0^2} ) = 200$ (red) and 
$-200$ (blue) in time-averaged flow field.
Figure 7(c) demonstrates that, similar to the case with resolution 
$\Delta r ^* = 23$, various modes of cavity surface oscillations are present along the entire length of the cavity in simulations when mesh resolution 
$\Delta r^* = 18.5$ is employed. This comparison indicates that the coarsest mesh is unable to capture the cavity surface oscillations, and these oscillations are present in the results of both cases with 
$\Delta r ^* = 23$ and 
$18.5$. These results demonstrate that in addition to the necessity of an adequately fine mesh for accurate simulation of the non-cavitating tip vortex flow, such mesh resolutions are crucial for the simulations to capture the cavity surface oscillations of tip vortex.
4.2.2. Overview of tip vortex cavity surface oscillations
The cavity surface oscillations are examined by employing a method similar to what is commonly used in experimental studies (Pennings et al. Reference Pennings, Bosschers, Westerweel and van Terwisga2015a) to extract various modes of cavity surface oscillations, which uses the diameter of the cavity viewed from different sides to distinguish between different modes of oscillation. Furthermore, the results obtained using this method are compared with the analytical solution discussed in Bosschers (Reference Bosschers2018b).
In order to extract various modes of the cavity surface oscillations, illustrated in figure 1, from the complex cavity surface similar to the tip vortex cavity illustrated in figure 7, the tip vortex cavity diameter observed from the top and side views is employed in the experimental study of Pennings et al. (Reference Pennings, Bosschers, Westerweel and van Terwisga2015a). This method can be utilized in numerical studies to attain an overview of the oscillations captured in the simulations, which can be used as a basis for more complex post-processing of the results. The non-dimensional tip vortex cavity radius observed from the side (
$xz$ view) and top (
$xy$ view) of the domain is presented in figure 9, which clearly indicates the presence of cavity surface oscillations in the results of the simulations. The radius used for non-dimensionalizing the data throughout this work, i.e. 
$r_c$, is the average of the spatially and temporally averaged radius observed from the top and side views in the region of interest (
$0.1 < x/c_0 < 1.2$).
Figure 9. Spatial-temporal variation of the non-dimensional radius observed at 
$\sigma = 1.7$ from (a) side view at the 
$xz$-plane, and (b) top view at the 
$xy$-plane.
Figure 9 shows that the cavity surface exhibits a relatively stationary behaviour from the tip until 
$x/c_0 \sim 0.5$ (formation region), and transitions to a more dynamic behaviour further downstream (developed region). Comparing figures 9(a) and 9(b) reveals that the spatial diameter variations of the roll-up region of the cavity (
$0.1 < x/c_0 < 0.5$) observed from these two views are out of phase, which indicates the presence of a nearly stationary twisting shape (double-helical mode) in this region of the cavity, similar to the observations of Ye et al. (Reference Ye, Wang and Shao2023). This stationary double-helical deformation is also evident in the time-averaged tip vortex cavity illustrated in figure 8(a). In order to extract the cavity oscillation modes from the side-view and top-view diameters, similar to the method employed by Bosschers (Reference Bosschers2018b), the cross-power spectral density (CPSD) is calculated from the two-dimensional fast Fourier transform (FFT) of the data obtained from the two views using the equation
\begin{equation} S(k_x, f) = 120 + 10 \log_{10} \left[ \frac{G_{top}(k_x, f)\,G^*_{side}(k_x, f)}{r_c^2} \right] , \end{equation}
where 
$k_x$ denotes the streamwise wavenumber, 
$G(k_x, f)$ represents the 2-D FFT of the radius data observed from the corresponding view, and 
$G^*(k_x, f)$ is its complex conjugate.
The wavenumber–frequency dependence of 
$S(k_x, f)$ is depicted in figure 10(a), where the frequency and wavenumber are non-dimensionalized using the flow characteristic frequency scale 
$f_{ref} = {U_\infty }/{c_0}$ and the root chord length 
$c_0$, respectively. Various modes of cavity oscillation manifest themselves as curves with higher 
$S(k_x, f)$ values, as shown in figure 10(a). The phase difference between the 2-D FFT of the radius data obtained from the side and top views is also required in order to distinguish between the breathing and double-helical modes. The phase difference, illustrated in figure 10(b), is nearly equal to 
$180^{\circ }$ in the double-helical mode, and 
$0^{\circ }$ in breathing mode oscillations. The results depicted in these figures clearly indicate the presence of the breathing and the double-helical modes.
Figure 10. Cavity oscillatory mode extraction using two orthogonal views of the cavity. (a) The CPSD obtained from (4.1). (b) Phase difference between 2-D FFT of cavity diameter viewed from the side and top.
The semi-analytical solution of Bosschers (Reference Bosschers2018b) is also plotted in figure 10(a) for comparison with the results obtained from the simulations in this work. The analytical solution for the dispersion relation of the cavity surface waves is
\begin{equation} \frac{2 {\rm \pi}r_c f ^ {{\pm}}}{U_\infty} = \tilde{U}_x k_x r_c + \tilde{U}_\theta n \pm \sqrt{K_\sigma} \sqrt{\frac{- | k_x r_c | \, K'_n(| k_x r_c |)}{K_n(| k_x r_c |)}}\,T_\omega , \end{equation}
where 
$\tilde {U}_x$ and 
$\tilde {U}_\theta$ are the non-dimensional axial and azimuthal velocities, respectively, 
$K_n$ and 
$K'_n$ denote the modified Bessel function of the second kind and its first derivative, respectively, 
$K_\sigma$ is a non-dimensional stiffness coefficient, and 
$T_\omega$ represents the effect of surface tension, which is neglected in this work (
$T_\omega = 1$).
According to (4.2), accurate prediction of the velocity profile of the tip vortex flow under cavitating conditions is vital for the cavity surface waves to be captured correctly by numerical simulations, which indicates the necessity of comparison of the velocity profiles in cavitating conditions with experimental data. The azimuthal velocity profiles obtained from simulations with different cavitation numbers are compared with the experimental results reported by Pennings et al. (Reference Pennings, Westerweel and Terwisga2015b) in figure 11. A maximum error of 
$2.2\,\%$ is observed between the maximum non-dimensional velocity values obtained from the simulations and the experimental data considering the uncertainty of the experimental results, which demonstrates that the azimuthal velocity profiles obtained from the simulations conducted in this work agree well with the experimental data for various cavitation numbers.
Figure 11. Comparison of non-dimensional azimuthal velocity profiles in simulations with experimental results of Pennings et al. (Reference Pennings, Westerweel and Terwisga2015b).
It is evident from the experimental results of Pennings et al. (Reference Pennings, Westerweel and Terwisga2015b) that the occurrence of cavitation influences the azimuthal velocity profile of the wetted tip vortex by decreasing the maximum azimuthal velocity and shifting its location towards larger radii. Bosschers (Reference Bosschers2018a) proposed a model for the azimuthal velocity profile of cavitating tip vortices as
\begin{equation} U_\theta (r) = \frac{\varGamma}{2 {\rm \pi}r} \left[ 1 - \frac{r_v^2}{r_v^2 + \zeta r_c^2} \exp \left( -\zeta\,\frac{r^2 - r_c^2}{r_v^2} \right) \right] , \end{equation}
where 
$\zeta = 1.2564$. Based on (4.3), the azimuthal velocity profile exhibits a local maximum where
\begin{equation} \frac{\partial U_\theta}{\partial r} = \frac{\varGamma}{2 {\rm \pi}r^2} \left[{-}1 + \frac{r_v^2}{r_v^2 + \zeta r_c^2} \left( 1 + 2 \zeta\,\frac{r^2}{r_v^2} \right) \exp \left( -\zeta\,\frac{r^2 - r_c^2}{r_v^2} \right) \right] = 0 , \end{equation}
which indicates that the radial position of the maximum azimuthal velocity in a cavitating tip vortex depends merely on the cavity radius 
$r_c$, which represents the extent of cavitation occurring within the tip vortex, and the wetted tip vortex radius 
$r_v$. On the other hand, it is evident in figure 11 that the radial position of the maximum azimuthal velocity in the cavitating tip vortex flow field in the simulations agrees with the experiments. It is concluded that the extent of cavitation within the tip vortex flow is accurately predicted in the simulations.
The coefficients in the analytical solution, i.e. 
$\tilde {U}_x$, 
$\tilde {U}_\theta$ and 
$K_\sigma$, are determined using the results of the simulations and by fitting the analytical solution onto the dispersion relation curves obtained from the simulations. The calculation of non-dimensional velocities 
$\tilde {U}_x$ and 
$\tilde {U}_\theta$ from the numerical results is carried out through spatial averaging of the axial and azimuthal velocity values on the iso-surface of 
$p=p_{sat}$ in the tip vortex flow field, since this is used as a boundary condition to derive the analytical solution. As confirmed by Bosschers (Reference Bosschers2018b), the stiffness coefficient 
$K_\sigma$ can be approximated as 
$\tilde {U}_\theta ^ 2$. This method yields 
$\tilde {U}_x = 1.03$, 
$\tilde {U}_\theta = 0.61$ and 
$K_\sigma = 0.37$ for the coefficients. Tuning these coefficients to fit the analytical solution to the numerical dispersion relations yields 
$\tilde {U}_x = 1.05$, 
$\tilde {U}_\theta = 0.68$ and 
$K_\sigma = 0.32$, which agree well with the values calculated using the numerical results of the velocity profiles.
The results illustrated in figure 10(a) indicate that all of the four double-helical mode oscillations (
$n=\pm 2^{\pm }$) are captured in the simulations, and these modes of cavity oscillation observed in the simulations agree reasonably well with the analytical results. The modes 
$n = 2^+$ and 
$n = -2^{-}$ observed in the results obtained from the simulations deviate from the analytical results, which might be resolved if the parameters of (4.2) are adjusted further. Moreover, some high-CPSD curves with 
$0^{\circ }$ phase difference can be observed in figure 10(a) that are not predicted by the analytical model.
Another dominant stripe observed in figure 10(a), which exhibits a phase difference of 
$0^{\circ }$ between the two views, shows the convection of the disturbances on the cavity surface. The convection group velocity (
$c_g$) in this case is found to be 
$c_g = 0.9 U_\infty$, which is similar to the value used by Bosschers (Reference Bosschers2018b), i.e. 
$c_g = 0.95 U_\infty$. Moreover, the dispersion curve corresponding to the 
$n=1^-$ mode collapses onto the convection line, as observed in figure 10(a). In addition, figure 10(a) shows that the wavenumber–frequency lines obtained for the breathing and double-helical modes are close in this case, and the dominance of the double-helical mode oscillations in this case hinders the detection of the breathing mode oscillations. Only traces of the breathing mode oscillations can be seen in the phase difference between the two views shown in figure 10(b), which demonstrates that a more effective method is required for the extraction of the breathing mode from the cavity oscillations. These results indicate that the numerical framework developed in this work is capable of accurately capturing the tip vortex cavity surface oscillations, which allows further investigation of the cavity surface oscillations and their contribution to the radiated noise.
4.3. Breathing mode oscillations
A major proportion of the TVC noise, especially at low and medium frequencies, is attributed to the mode of oscillation responsible for cavity volume variations, i.e. the breathing mode (Wang et al. Reference Wang, Bai, Cheng, Ji and Peng2023). Moreover, the centre frequency of the broadband hump observed in the URN spectrum is hypothesized to be related to the resonance frequency of the tip vortex cavity (Bosschers Reference Bosschers2018b), and experimental evidence has been provided by Pennings et al. (Reference Pennings, Bosschers, Westerweel and van Terwisga2015a) showing that this resonance frequency is the frequency of the zero group velocity of the breathing mode. These connections linking the breathing mode to the TVC noise necessitate a comprehensive investigation of this mode of cavity oscillation. Therefore, in this subsection, this mode of oscillation is investigated in detail for three cases with different cavitation numbers. First, in § 4.3.1, a method is developed to extract the breathing mode from the cavity surface obtained from the simulations by defining an effective radius for the tip vortex cavity. The temporally and spatially averaged effective radii of the cavity are analysed in §§ 4.3.2 and 4.3.3, respectively. Subsubsection 4.3.4 focuses on the spatial and temporal development of the effective radius, delving into the characteristics of local breathing mode oscillations. The POD is performed on the effective radius data in § 4.3.5 to further analyse the streamwise development of the tip vortex cavity's volume.
4.3.1. Breathing mode extraction
In order to examine the breathing mode oscillations in detail, a method is required for extracting this mode of cavity oscillation. Klapwijk et al. (Reference Klapwijk, Lloyd, Vaz, van den Boogaard and van Terwisga2022) fitted an ellipse onto the transverse cross-section of the cavity interface, and defined an effective radius as 
$r_{eff}=\sqrt {ab}$, where 
$a$ and 
$b$ are the semi-major and semi-minor axes of the fitted ellipse, respectively. However, as explained by Wang et al. (Reference Wang, Bai, Cheng, Ji and Peng2023), the method used in Klapwijk et al. (Reference Klapwijk, Lloyd, Vaz, van den Boogaard and van Terwisga2022) might lead to loss of deformation information when the cross-section is non-elliptical. Such a highly non-elliptical transverse cross-section of the cavity was also observed in our simulations, as shown in figure 12.
Figure 12. Calculation of the effective radius using the tip vortex cavity surface data.
Since the breathing mode is the only cavity oscillation mode responsible for the cavity volume variations, the idea of employing the cavity's cross-sectional surface area for calculating an effective radius is a sound approach. Based on this approach, first, the transverse cross-section of the cavity at different downstream locations is extracted, the area enclosed by which is then calculated. Subsequently, an effective radius is calculated based on the area enclosed by the cavity cross-section as 
$A_{enclosure} = {\rm \pi}r_{eff}^2$, as depicted in figure 12. The effective radius obtained using this method contains only the volume variation information, and other forms of oscillation are excluded from this quantity. Therefore, this parameter can be employed to investigate breathing mode oscillation in the absence of the interference of other modes.
The two-dimensional Fourier transform of the effective radius obtained using this method for the case with cavitation number 
$\sigma = 1.7$ is illustrated in figure 13, which displays the wavenumber–frequency dependency of the breathing mode cavity oscillations, together with the analytical solution. Based on this figure, the breathing mode oscillations obtained from performing the extraction method on the simulation data agree well with the analytical solution, which shows that this extraction method, together with the numerical framework employed in this work, is able to accurately capture the breathing mode oscillations.
Figure 13. Frequency–wavenumber diagram of the effective radius data at 
$\sigma = 1.7$.
4.3.2. Temporally averaged effective radius
The tip vortex cavity surface can be decomposed into a mean cavity surface and time-dependent cavity surface oscillations. The formation of the mean cavity interface is related to the overall roll-up process of the tip vortex. In this subsubsection, the impact of the roll-up process on the mean cavity volume is investigated by analysing the time-averaged effective radius. The goal here is to gain insight into the evolution of the tip vortex cavity from its inception to decay.
Figure 14 shows the temporally averaged effective radius, denoted by 
$\langle r_{eff}\rangle$ and non-dimensionalized by 
$r_c$, with respect to the streamwise location (
$x$) non-dimensionalized by the root chord length (
$c_0$) for three cavitation numbers. According to the streamwise development of the temporally averaged effective radius in these cases, the cavity undergoes a region of growth, where the cavity is forming and the cavity effective radius has a growing trend overall, followed by a decay region.
Figure 14. Streamwise evolution of the temporally averaged effective radius at different cavitation numbers.
A comparison of different cavitation numbers indicates that the growth region shrinks as the cavitation number increases. In the case with 
$\sigma = 2.6$, the decay begins approximately at 
$x/c_0 = 0.12$. In the case with 
$\sigma = 1.7$, the growth region extends until 
$x/c_0 = 0.4$, and when 
$\sigma = 1.2$, the growth continues further until 
$x/c_0 = 0.53$. Furthermore, figure 14 reveals that decay occurs at a greater slope in higher cavitation numbers. The most sustained cavity is observed in the case with 
$\sigma = 1.2$, where the slope of decay (
${{\rm d}r^*}/{{{\rm d}{\kern0.8pt}x}^*}$, where 
$r^* = {r_{eff}}/{r_c}$ and 
$x^* = {x}/{c_0}$) is approximately 
$0.146$. This decay slope in the cases 
$\sigma = 1.7$ and 
$2.6$ increases to 
$0.288$ and 
$1.027$, respectively.
An interesting feature of the temporally averaged effective radius seen in figure 14 is the spatial oscillations observed in the growth region. In the case with 
$\sigma = 1.7$, the cavity volume exhibits periodic spatial variation with the downstream location in the growth region (
$x/c_0 < 0.4$); however, in the decay region, these spatial variations of the cavity volume disappear. When 
$\sigma = 2.6$, such variations are present until 
$x/c_0 = 0.23$, which can be attributed to the smaller growth region in this case compared to 
$\sigma = 1.7$. On the other hand, as is evident in this figure, the spatial volumetric variations in the formation region vanish when the cavitation number is lowered to 
$\sigma = 1.2$.
The spatial variations of the temporally averaged effective radius in the cavity growth region are further investigated by calculating their Fourier transform, which is illustrated in figure 15. The horizontal axis in this figure represents the streamwise wavenumber (
$k_x$) non-dimensionalized by the root chord length (
$c_0$). The large-amplitude low-wavenumber components (
$k_x c_0 < 10$) correspond to the overall growth and decay of the cavity. An interesting feature observed in the case 
$\sigma = 1.7$ is the hump within the dimensionless wavenumber range 
$80 < k_x c_0 < 177$. The centre wavenumber of this hump is 
$k_x c_0 = 100$, which corresponds to wavelength 
$\lambda = 0.0628 \, c_0$. A local maximum is also evident for the case 
$\sigma = 2.6$, which is less significant than the hump observed in the case 
$\sigma = 1.7$ due to the spatial oscillations decaying earlier in the 
$\sigma = 2.6$ case compared to 
$\sigma = 1.7$. The wavenumber of this hump is similar to that of the hump observed in the 
$\sigma = 1.7$ case, indicating that these spatial variations originate from the tip vortex flow features and are not related to the cavity characteristics. Moreover, the existence of this component in both cases further supports the idea that it stems from the characteristics of the tip vortex flow rather than the cavity's characteristics.
Figure 15. Fourier transform of the non-dimensional temporally averaged effective radius variations.
To further investigate this phenomenon, the time-averaged flow field is analysed at cavitation number 
$\sigma = 1.7$. The time-averaged tip vortex cavity depicted in figure 8(a) exhibits a stationary double-helical structure in the growth region. The iso-surfaces of the streamwise component of vortex stretching, i.e. 
$(\boldsymbol {\omega }\boldsymbol {\cdot }\boldsymbol {\nabla }) u$, normalized by 
${U_\infty ^2}/{c_0^2}$ and illustrated in figure 8(b), exhibit twisted patterns similar to the time-averaged cavity structure. This similarity demonstrates that the roll-up of streamwise coherent structures from the boundary layer into the tip vortex generates negative and positive vortex stretching regions swirling around the tip vortex cavity in a helical pattern, which leads to the stationary twisted shape of the tip vortex cavity within the growth region. The interactions between these structures and the tip vortex cavity may be a source of the spatial volume variations observed within the growth region of the tip vortex cavity.
The possible range of pitch values for such helical structures is calculated to further examine the correlation between the coherent structures swirling around the tip vortex cavity and the cavity's stationary morphology. Such a coherent structure rolling up into the main tip vortex at radial distance 
$r$ from the cavity centre would have pitch value 
$({2 {\rm \pi}r}/{U_\theta } )U_x$. The results presented in figure 11 indicate that such a coherent structure would have a streamwise-to-azimuthal velocity ratio within the approximate range 
$1 < {U_x}/{U_\theta } < 2$, located within 
$r_v$ to 
$2r_v$ from the centre of the cavity. It can be concluded that the pitch of such helical coherent structures would be within a range that corresponds to the non-dimensional streamwise wavenumber range 
$k_x c_0 \in [29, 114]$. It is evident in figure 15 that the high-amplitude components of the spatial cavity volume variations occur within the range 
$k_x c_0 \in [30, 180]$. Comparison of these two ranges demonstrates that the range of possible pitch values of coherent structures swirling around the tip vortex cavity and the range of wavenumbers at which the spatial cavity volume variations occur are in good agreement, further supporting the correlation between the coherent structures swirling around the cavity and the stationary structure of the cavity.
4.3.3. Spatially averaged effective radius
The volume of the tip vortex cavity is an important quantity that can reveal any cyclic process that the cavity might undergo. In this subsubsection, the spatially averaged effective radius, shown as 
$\overline {r_{eff}}$, is used to analyse the volume variations of the entire vapour cavity in the region of interest (
$0.1 < x/c_0 < 1.2$). The spatially averaged effective radius non-dimensionalized using 
$r_{c}$ is depicted in figure 16(a) as a function of the non-dimensional time (
$t^* = {t U_\infty }/{c_0}$) for the three cavitation numbers. The spatial averaging for the case 
$\sigma = 2.6$ is carried out for the streamwise location range 
$0.1 < x/c_0 < 1.0$ because the cavity collapses between 
$x/c_0 = 1.0$ and 
$x/c_0 = 1.2$ in this case. It is evident in figure 16(a) that the volume of the cavity in the region of interest undergoes temporal variations, which occur more intensely when 
$\sigma = 1.7$, where cavitation occurs moderately, compared to the other cases.
Figure 16. Temporal dynamics of spatially averaged effective radius. (a) Non-dimensional spatially averaged effective radius temporal development. (b) Spectrum of the non-dimensional spatially averaged effective radius variations.
The frequency spectra of the spatially averaged effective radius fluctuations are illustrated in figure 16(b) obtained using the FFT algorithm. The spectra of the fluctuations demonstrate that a large proportion of the vapour volume variations occur in the low-frequency region 
${f c_0}/{U_\infty } < 3$, which approximately corresponds to the range of frequencies below 
$165 \ \mathrm {Hz}$. Moreover, these spectra exhibit specific features that require further examination. The dimensionless frequency 
${f c_0}/{U_\infty }$ is denoted by 
$f^*$ in the rest of this paper.
In the case 
$\sigma = 1.2$, a local maximum is observed at 
$f^* = 0.37$, and there are peaks observed in the spectrum corresponding to the 3rd, 4th and 6th harmonics of this frequency, i.e. 
$f^* = 1.11$, 1.48 and 2.22, respectively, with the 3rd harmonic at 
$f^* = 1.11$ being the most significant. Furthermore, the spectrum exhibits a local maximum at the subharmonic of this distinct peak at 
$f^* = 0.55 \sim {1.11}/{2}$. Since there is no other significant fluctuating component evident in the spectrum in this case, it can be concluded that all of the major components of this spectrum originate from the same underlying source.
In the spectrum obtained from the case 
$\sigma = 1.7$, a peak is evident at 
$f^* = 0.42$, together with its 2nd, 4th and 5th harmonics, corresponding to 
$f^* = 0.83$, 1.66 and 2.12, respectively. The spectrum also exhibits a significant fluctuating component at 
${f c_0}/{U_\infty } = 0.23$, which is not shown in the spectrum since its corresponding dimensional frequency falls below 20 Hz. Furthermore, a distinct peak is observed at 
$f^* = 1.02$ in this case. In addition, the spectrum in this case displays a peak at 
$f^* = 0.55$, similar to the 
$\sigma = 1.2$ case, and its 4th harmonics at 
$f^* = 2.22$. Another peak is also observed in this case, at 
$f^* = 1.15$, which can be the 2nd harmonic of the 
$f^* = 0.55$ component.
In the case with cavitation number 
$\sigma = 2.6$, the time-averaged effective radius does not experience significant fluctuations, with local maxima at 
$f^* = 0.46$ and its 2nd harmonic 
$f^* = 0.92$. Furthermore, the component at 
$f^* = 2.22$ is observed in this case in a similar way to the other cavitation numbers, in addition to its 2nd harmonic at 
$f^* = 4.43$. In this case, a peak is observed at 
$f^* = 4.25$ that is not evident in other cases.
The features of the spatially averaged effective radius spectra discussed herein indicate common features among different cavitation numbers. The presence of such components for all cavitation numbers indicates that these fluctuation components originate from the flow characteristics. In addition, in all of the cases, it is evident that the most intense oscillations of the spatially averaged effective radius occur at frequencies close to the flow characteristic frequency. On the other hand, the spectra of spatially averaged effective radius fluctuations of these cases differ in some other features. Such different components among different cavitation numbers indicate that these components might be due to the specific characteristics of the cavity that forms at each cavitation number.
Since the strength of the tip vortex is dependent on the lift coefficient, the lift fluctuations potentially affect the overall variations of the tip vortex cavity volume. Therefore, the temporal variations of the lift coefficient are investigated to reveal their potential impact on cavity volume variations. The lift coefficient as a function of time obtained from the simulations is depicted in figure 17(a), and figure 17(b) illustrates its frequency spectrum. As is evident, the spectrum of the lift coefficient exhibits a distinct peak at 
$f^* = 0.21$, and less significant local maxima at 
$f^* = 0.35$ and 
$0.48$.
Figure 17. Temporal dynamics of lift coefficient. (a) Temporal variations of lift coefficient. (b) Frequency spectrum of lift coefficient variations.
Comparison of the fluctuating components of the lift coefficient with the features of spatially averaged effective radius spectra reveals some similarities between them. In the case 
$\sigma = 1.7$, the most distinct peak of spatially averaged effective radius fluctuations occurred at frequency 
$f^* = 0.23$, which is similar to the major fluctuating component of the lift coefficient spectrum. Furthermore, the case 
$\sigma = 1.2$ exhibited mainly the harmonics of 
$f^* = 0.37$ in its volume variation spectrum, which is consistent with the second peak observed in the lift coefficient spectrum. The peak observed at 
$f^* = 0.42$ in the 
${\overline {r_{eff}}}/{r_c}$ spectrum of 
$\sigma = 1.7$ could also be related to this component of the lift coefficient fluctuations.
The results presented in this subsubsection show that the variations in the vapour volume of the tip vortex cavity are dominated mainly by the flow characteristics, and the characteristics of the cavity forming in different cases have less significant impact on the overall variations in the vapour volume.
4.3.4. Spatial-temporal characteristics of breathing mode oscillations
The effective radius calculated using the method explained in § 4.3.1 is a function of time and streamwise location; therefore, it can reveal the spatial and temporal evolution of the breathing mode oscillations. The spatial and temporal variations of the non-dimensionalized effective radius for the cavitation number 
$\sigma = 1.7$ are illustrated in figure 18(a), and figure 18(b) displays the non-dimensional fluctuations of the effective radius at every streamwise location by subtracting the time-averaged effective radius (
$\langle r_{eff}\rangle$) at the corresponding location from the effective radius value. The spatial variations observed in this case in the cavity growth region 
$0.1 < x/c_0 < 0.4$ can be observed in figure 18(a).
Figure 18. Spatial-temporal variations of the non-dimensional effective radius and its fluctuations (
$\sigma = 1.7$): (a) 
${r_{eff}}/{r_{c}}$, (b) 
$({r_{eff} - \langle r_{eff}\rangle })/{r_{c}}$.
Furthermore, it is evident in figure 18(b) that the behaviour of the effective radius fluctuations alters at 
$x/c_0 \sim 0.3$ in this case. In the region upstream of 
$x/c_0 = 0.3$, pairs of out-of-phase oscillating regions are observed marked by red dashed ellipses in figure 18(b). These pairs are potential evidence of adjacent antinodes of standing breathing mode waves, more evidence for which is provided in figure 19. Moreover, it can be observed in figure 18(a) that at streamwise location 
$x/c_0 \sim 0.22$, the effective radius periodically peaks approximately every 
$t_{ref} = {c_0}/{U_\infty }$ (flow characteristic time scale), indicating a fluctuating component of the effective radius at frequency 
$f_{ref} = {U_\infty }/{c_0}$ ( 
$f^* = 1$) at this location. On the other hand, in the region downstream of 
$x/c_0 =0.3$, the change in the oscillatory behaviour is apparent, where the propagation of disturbances is characterized by downstream propagation, as shown by a dashed line in figure 18(b). These lines, with a slope nearly equal to 1, display the convection of disturbances at a propagation speed close to 
$U_\infty$. Overall, these different observations within different regions of the cavity reveal the change in cavity oscillatory behaviour as it advances from the growth region to the decay region.
Figure 19. Spectrum of the effective radius fluctuations with respect to downstream location: (a) 
$\sigma =1.2$, (b) 
$\sigma =1.7$, (c) 
$\sigma =2.6$.
The spectra of the effective radius fluctuations as functions of streamwise location at different cavitation numbers are shown in figure 19. It is evident in this figure that similar to the behaviour of the spatially averaged effective radius fluctuations presented in § 4.3.3, the major proportion of the local fluctuations of the effective radius occurs at low frequencies (below 
$f^* = 5$). Furthermore, this figure further substantiates that the oscillation characteristics of the breathing mode alter as the cavity travels away from the tip. The presence of adjacent high- and low-amplitude components within the growth region in all of the cases, which are specified using magenta dashed ellipses, further supports the occurrence of standing waves within this region since these adjacent high- and low-amplitude regions represent the antinodes and nodes of standing breathing mode waves, respectively. Overall, according to figure 19, the breathing mode oscillations occur at lower frequencies in the growth region, and higher frequency components of the breathing mode oscillation intensify as the cavity travels downstream. The frequency spectra of breathing mode oscillations at three streamwise locations at different cavitation numbers are illustrated in figure 20. In the case 
$\sigma = 1.7$, the spectrum of the breathing mode oscillations at 
$x/c_0 = 0.1$ is also plotted, which reveals that the breathing mode oscillations are small in amplitude at this location, and the intensification of the breathing mode oscillations is apparent when the spectrum at this location is compared to locations further downstream. This indicates that these breathing mode oscillations are not initially present when the cavity forms, and are introduced onto the cavity surface as it flows downstream.
Figure 20. Frequency spectrum of the effective radius oscillations at multiple streamwise locations: (a) 
$\sigma =1.2$, (b) 
$\sigma =1.7$, (c) 
$\sigma =2.6$
In the case 
$\sigma = 1.7$, as mentioned earlier, adjacent regions of high- and low-amplitude fluctuations are observed within the growth region, with the main peaks occurring at 
$f^* = 1.02$, 1.15, 1.06 and 0.83. Within the transition region between the growth and decay regions, other fluctuating components, especially at higher frequencies, appear in the spectrum, such as the component at 
$f^* = 2.77$ and its 2nd harmonic 
$f^* = 5.54$. As the cavity reaches the decay region, the effective radius fluctuations exhibit a more broadband behaviour, with high-amplitude components at 
$f^* = 0.23$ and its 5th harmonic 
$f^* = 1.15$. Comparison of these results with the results presented in § 4.3.3 indicates that the fluctuations of the spatially averaged effective radius, which represents the volume of the cavity, are dominated by the local breathing mode oscillations within the growth region of the cavity.
In the case 
$\sigma = 1.2$, the transition to more broadband effective radius fluctuations is more pronounced when the cavity travels from the growth region to the decay region. The growth region of the cavity in this case displays spatially periodic high- and low-amplitude regions similar to those in other cases. These oscillations occur more intensely at 
$f^* = 0.55$ and its 2nd harmonic 
$f^* = 1.1$. Within the transition region, i.e. 
$0.4 < x/c_0 < 0.6$, other components in 
$f^* = 1.66$ and 
$f^* = 3.05$ are introduced into the breathing mode spectrum. In this case, the peak at 
$f^* = 0.55$ is present along the entire length of the cavity.
At cavitation number 
$\sigma = 2.6$, the oscillatory behaviour alters from that of lower cavitation numbers. The adjacent high- and low-amplitude regions within the growth regions at very low frequencies, 
$f^* < 0.2$, display a streamwise variation similar to that of 
$\sigma = 1.7$ within this region and frequency range. Furthermore, adjacent high- and low-amplitude regions are observed at frequency 
$f^* = 0.46$ within the growth region as well, although they also exhibit lower intensity. However, no significant broadband fluctuation is observed along the cavity, except for the end region, 
$x/c_0 > 0.9$. An interesting feature of the cavity in this case is the distinct fluctuating components within the location range 
$0.3 < x/c_0 < 0.5$ at frequencies within the range 
$4.53 < f^* < 4.71$.
Since the fine mesh region encompasses the entire cavity in the case 
$\sigma = 2.6$, the end effects of the tip vortex cavity can be examined in this case. The results indicate that the cavity undergoes a cyclic pattern of growth and detachment at its front, as illustrated in figure 21. The cycle starts with the elongation of the cavity, which progresses until it reaches a critical length, followed by detachment of a part of the cavity and shedding downstream. After detachment, the end of the tip vortex cavity retracts until it reaches a minimum length, and the cycle starts over with elongation. The frequency of this cycle was found to be within the range 
$4.5 < {f c_0}/{U_\infty } < 5.5$, which is similar to the peak frequencies observed in the breathing mode oscillation spectrum in figure 19(c) between 
$x/c_0 = 0.3$ and 
$x/c_0 = 0.6$. This highlights a potential correlation between the breathing mode oscillation within this frequency range and the cyclic pattern of growth and detachment observed in the tip vortex cavity's front.
Figure 21. The cyclic behaviour of the tip vortex cavity's front at 
$\sigma = 2.6$.
To further support that this cyclic tail behaviour is potentially related to the breathing mode of cavity surface oscillation rather than surface tension effects, an estimate of the Weber number is carried out using the equation
\begin{equation} We = \frac{\rho \overline{U_x}^2 r_c}{\gamma} , \end{equation}
where 
$r_c$ is the cavity radius close to the location of cavity breakup (
$x/c_0 = 1.0$), 
$\overline {U_x}$ is the azimuthally averaged axial velocity on the cavity surface at 
$x/c_0 = 1.0$, and 
$\gamma$ is the surface tension. Using (4.5), 
$We$ is estimated to be within the range 
$We \in [200, 300]$, which indicates that the inertial effects are dominant, and the effect of surface tension can be assumed to be negligible.
4.3.5. Decomposition of breathing mode oscillations
The analysis of the effective radius reveals that the breathing mode oscillations exhibit various characteristics in different regions of the cavity. In this subsubsection, POD is employed to decompose the effective radius data and obtain the main breathing mode oscillations. The time-averaged effective radius is subtracted from the results to extract the fluctuating components. The analysis in this subsubsection is carried out only for the case 
$\sigma = 1.7$.
The energy contribution of the effective radius fluctuation modes is plotted in figure 22(a), which shows that the fluctuations are distributed over a large number of modes; however, since the goal here is to gain insight into the fluctuations and not reduced-order modelling, the first few modes can be employed. Figure 22(b) shows the first 10 modes of effective radius fluctuations. The change in behaviour downstream of the cavity growth region, which corresponds to 
$x > 0.4$, is apparent in the fluctuation modes. It is also evident from the modes that the fluctuations in the cavity growth region have relatively small wavelengths and amplitudes compared to those of the decay region.
Figure 22. The POD of breathing mode oscillations (
$\sigma = 1.7$). (a) Singular values of the effective radius fluctuation modes. (b) First 10 modes of effective radius fluctuations obtained using POD (where ‘M’ denotes mode).
The amplification and alteration in the behaviour of the effective radius fluctuations downstream of the cavity growth region indicate that some oscillations are introduced onto the cavity interface in the roll-up process, as mentioned before. The source of these perturbations should be sought in the temporal variations of the tip vortex flow field and the roll-up of the wake into the tip vortex flow, which is not within the scope of this study.
4.4. Effect of cavity surface oscillations on pressure fluctuations
In order to investigate the contribution of cavity surface oscillations to the pressure fluctuations in different cases, the pressure is probed at 19 points within the domain, as shown in figure 23, on which an FFT is performed subsequently to obtain the spectra of the pressure fluctuations. The spectra of the pressure fluctuations at these points exhibit a similar behaviour in all cases; therefore, the pressure fluctuations spectrum at only one point, which is located 
$0.5c_0$ away from the tip in the spanwise direction (
$x=0$, 
$y=0$, 
$z=0.5c_0$), is reported herein and is shown in figure 24 for the three cases.
Figure 23. Locations of pressure probing.
Figure 24. Spectrum of the pressure fluctuations at point 
$x=0$, 
$y=0$, 
$z=0.5c_0$.
It is evident in figure 24 that decreasing the cavitation number leads to more intense pressure fluctuations. This indicates that tip vortex cavitation is, indeed, contributing to the pressure fluctuations within the domain. Moreover, it can be observed that as the cavitation number decreases, the spectrum shifts towards lower frequencies overall. Fluctuations in pressure at low frequencies (below 
$f^* = 1$) are relatively significant in all cases. In addition to these general features, the pressure fluctuations spectrum of each case exhibits some unique characteristics.
The spectrum of pressure fluctuations for 
$\sigma = 1.2$ shows a distinct peak at 
$f^* = 2.22$, with less significant peaks at its subharmonics 
$f^* = 1.11$ and 
$0.55$. As discussed in § 4.3.3, these components were observed in the spatially averaged effective radius of the cavity in the 
$\sigma = 1.2$ case. In addition, the presence of these components in local breathing mode oscillations was also discussed in § 4.3.4. It is observed in figure 20(a) that the peaks at 
$f^* = 1.11$ and 0.55 are more significant within the growth region, and the higher frequency peak at 
$f^* = 2.22$ is distinct at 
$x/c_0 = 0.98$ within the decay region. Therefore, it can be concluded that the growth and decay regions of the cavity contribute differently to the pressure fluctuations, with the contribution of the decay region being more significant since the peak at 
$f^* = 2.22$ possesses the highest amplitude in the pressure fluctuations spectrum. Furthermore, the pressure fluctuations spectrum in this case exhibits other peaks at 
$f^* = 0.74$, which is the 2nd harmonic of 
$f^* = 0.37$ observed in the spatially averaged effective radius spectrum, and 
$f^* = 1.66$, the 3rd harmonic of 
$f^* = 0.55$.
At cavitation number 
$\sigma = 1.7$, the pressure fluctuations spectrum shows peaks at low frequencies such as 
$f^* = 0.42$, and less significantly at 
$f^* = 0.55$, 0.65, 0.79 and 0.92. The low-frequency components at 
$f^* = 0.42$, 0.55 and 0.69 were also observed in the spatially averaged effective radius in this case in figure 16(b). The component at 
$f^* = 0.42$ is also evident within the growth region in figure 19(b), and the component at 
$f^* = 0.79$ is observed within the transition region from the growth regions to the decay regions, i.e. 
$0.4 < x/c_0 < 0.6$. In the pressure fluctuations spectrum, in this case, a hump is observed within the frequency range 
$3.4 < f^* < 5.1$, which is the peak at 
$f^* = 3.79$. This peak is also observed in figure 19(b) within the decay region. Moreover, this peak is the 9th harmonic of 
$f^* = 0.42$ observed in the breathing mode spectrum within the growth region. Overall, similar to the 
$\sigma = 1.2$ case, the results demonstrate the different contributions of the breathing mode of oscillation within different regions of the cavity to the pressure fluctuations.
In the case 
$\sigma = 2.6$, the pressure fluctuations spectrum exhibits a relatively small hump within the frequency range 
$4.1 < f^* < 5.4$, with a peak at 
$f^* = 4.71$. The spectrum of breathing mode oscillations for this case depicted in figure 19(c) revealed intense breathing mode oscillations within the frequency range 
$4.53 < f^* < 4.71$ within the streamwise location range of 
$0.28 < x/c_0 < 0.55$. Moreover, the spatially averaged effective radius spectrum presented in figure 16(b) exhibited a peak at 
$f^* = 1.66$, the subharmonic of which is observed in figure 24 at 
$f^* = 0.83$. Furthermore, the cyclic growth and detachment behaviour of the cavity's end was found to be occurring within a non-dimensional frequency range 
$4.5 < f^* < 5.5$, which indicates that this cyclic tail effect is potentially influencing the pressure fluctuations within the domain as well.
An interesting common characteristic is observed in these results, which indicates that the peak frequencies of the hump in the pressure fluctuations spectra, in addition to being related to breathing mode oscillatory components within the decay region, are potentially related to a harmonic of the frequency at which the adjacent low- and high-amplitude regions within the growth region oscillate. In cases with 
$\sigma = 1.2$, 
$1.7$ and 
$2.6$, the peak frequency of the hump was found to be 
$f^* = 2.22$, 3.79 and 4.71, respectively. These components correspond to the 4th, 9th and 10th harmonics of the fluctuating components observed within their corresponding cavities’ growth regions, which were 
$f^* = 0.55$, 
$0.42$ and 
$0.46$. In addition, it is observed that at higher cavitation numbers, the peak of the pressure fluctuations corresponds to a higher harmonic of the breathing mode oscillatory component within the growth region.
Based on the correlations between the cavity breathing mode oscillations and the pressure fluctuations, it can be concluded that the growth region of the cavity is correlated with the low-frequency peaks observed in the pressure fluctuations spectrum, and the mid-range of this spectrum, where the hump in the pressure fluctuations spectrum occurs as well, is correlated with the breathing mode oscillations in the decay region of the tip vortex cavity.
5. Concluding remarks
In this study, tip vortex cavitating flow over a stationary elliptical NACA66(2)-415 hydrofoil was investigated numerically using a recently developed finite-element-based cavitation flow solver for high Reynolds numbers. To investigate the effect of mesh resolution on the simulation results and establish a method for estimation of the required mesh resolution for large eddy simulations (LES) of TVC based on the flow properties, a pressure-gradient-based length scale was developed employing the Rankine vortex model, the Kutta–Jukowski theorem, and the turbulent boundary layer thickness. This new length scale was then utilized for non-dimensionalization of the mesh resolution within the tip vortex flow region. The simulations carried out in this study for various non-dimensional mesh resolution values (
$\Delta r ^ *$) revealed that a 
$\Delta r ^ *$ value of 18.5 or lower should be utilized for the accurate simulation of tip vortex cavitating flows using LES. The effect of mesh resolution on cavity surface oscillations was also investigated, the results of which indicated that resolution 
$\Delta r ^ * = 18.5$ was able to capture these oscillations as well. The oscillatory dynamics of the tip vortex cavity was investigated in the results obtained from the simulations carried out with a mesh satisfying the proposed resolution requirements. The simulations were conducted for three different cases that differed in cavitation number. The cavity dynamics captured in the simulations was shown to agree well with the analytical solution. In addition, the azimuthal velocity profile of the cavitating tip vortex agreed with the experimental data, further validating the numerical framework employed in this work.
Since it is known that the breathing mode of oscillation is the most influential mode of the surface oscillations concerning the generated noise compared to other modes, this mode of oscillation was extracted from the numerical results using a parameter defined based on the tip vortex cavity cross-sectional area, referred to as the cavity effective radius (
$r_{eff}$). Investigation of the temporally averaged effective radius demonstrated that the tip vortex cavity experiences different regions as it progresses away from the tip. The results indicated that the growth region – i.e. the region close to the tip where the cavity is forming and experiences an overall increasing trend in volume – shrinks with increasing cavitation number. Furthermore, it was demonstrated that the decay occurs at an increased slope in scenarios with higher cavitation numbers. Another observation was the occurrence of spatial periodic cavity volume variations within the growth region. These variations exhibit similar dominant wavenumbers between cases with different cavitation numbers, which reveals that the dominant wavenumber of such oscillations is related to the mean tip vortex roll-up process, and not the characteristics of the cavity forming in each specific case. Further analysis of the vortex stretching term within the tip vortex flow region indicated the presence of negative and positive vortex stretching regions swirling around the cavity forming the stationary double-helical shape of the tip vortex cavity. Moreover, the possible pitch values for the coherent structures swirling around the tip vortex cavity were found to agree with the wavenumber range of the spatial volume variations of the cavity, further supporting the argument that the roll-up of coherent structures into the tip vortex is responsible for the spatial cavity volume variations.
Investigation of the spatially averaged effective radius indicated that the overall vapour volume undergoes temporal variations, with the most intense oscillations occurring in the moderately cavitating case, i.e. 
$\sigma = 1.7$. Furthermore, it was observed that the most significant cavity volume fluctuation components possessed frequencies similar to the flow characteristic frequency. The temporal variations of the lift coefficient were examined, which revealed potential correlations between the lift force fluctuations and the tip vortex cavity overall volume variations. These analyses indicated that the flow characteristics play a major role in the variations of overall cavity volume.
Subsequently, the local breathing mode oscillations were investigated. The results provided evidence for the potential occurrence of waves in standing breathing mode within the growth region of the tip vortex cavity. More evidence was provided for the alteration of the cavity's oscillatory behaviour between the growth and decay regions, employing POD. Moreover, the results indicated that the breathing mode oscillations are not inherent to the tip vortex cavity, and such oscillations appear and intensify as the cavity travels away from the tip. Therefore, the tip vortex cavity displays different characteristics in the growth and decay regions. The breathing mode oscillations within the decay region were observed to extend over a more broadband range of frequencies compared to the growth region.
Overall, the breathing mode oscillations within the growth region were observed to occur within the low-frequency range, and such oscillations within the decay region were found to possess higher frequencies. The simulations revealed that at a low cavitation number, similar to the case 
$\sigma = 1.2$ in this work, no significant peak was observed in the breathing mode spectrum; however, distinct peaks were observed at higher cavitation numbers, with the most intense ones occurring in the case 
$\sigma = 1.7$. Furthermore, similar to the spatially averaged effective radius fluctuation spectrum, the local breathing mode oscillation in the case 
$\sigma = 2.6$ extended to higher frequencies compared to lower cavitation numbers. An interesting feature of the tip vortex cavity was observed when 
$\sigma = 2.6$, which was the cyclic growth and detachment behaviour of the cavity at its end. It was observed that the cavity undergoes a process of elongation, detachment and retraction at its end, which was found to occur within a frequency range 
$4.5 < {f c_0}/{U_\infty } < 5.5$. The local breathing mode oscillations in this working condition also showed fluctuations within this frequency range, indicating a potential correlation between the cyclic tail behaviour of the tip vortex cavity and its breathing mode oscillations. This aspect of tip vortex cavitation requires further investigation.
The pressure fluctuations within the domain were investigated in cases with different cavitation numbers. The results indicated that decreasing the cavitation number, which leads to a higher extent of TVC, results in more intense pressure fluctuations. Another effect of decreasing the cavitation number was found to be shifting the pressure fluctuations spectrum towards lower frequencies, which is similar to the behaviour of breathing mode oscillations. Comparison of the pressure fluctuations spectra at different cavitation numbers with the overall vapour volume variations and local breathing mode oscillations at the corresponding cavitation numbers demonstrated the existence of correlations between the pressure fluctuations and the breathing mode of oscillation. The low-frequency components of the pressure fluctuations were found to be correlated with the growth region of the cavity. The centre frequencies of the humps observed in the pressure fluctuations spectra were found to be related to the breathing mode oscillations within the decay region in all cases. Despite the existence of these correlations, there were some features in the pressure fluctuations that were not related to the breathing mode of the cavity surface oscillation. Moreover, some significant fluctuating components were observed in the breathing mode oscillations that did not affect the pressure fluctuations despite having high amplitudes. Therefore, further investigation is required to thoroughly understand the correlations between the breathing mode oscillations and pressure fluctuations.
Funding
This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through the grant NSERC IRCPJ 550069-19. The computational resources were provided by Compute Canada (www.alliancecan.ca) and the Advanced Research Computing facility at the University of British Columbia (www.arc.ubc.ca). We would like to thank Dr T. van Terwisga for providing the geometry file of the elliptical NACA66(2)-415 hydrofoil. S.L. thanks Mr N.B. Darbhamulla for useful discussions during the development of this work.
Declaration of interests
The authors report no conflict of interest.
