Finite Volume, Computational Fluid Dynamics-Based Investigation of Supercavity Pulsations

Due to their potential and proven hydrodynamic applications to high-speed drag reduction and facilitation of high-speed marine lifting surfaces, ventilated gas supercavities are of particular interest. The reduced drag is achieved when a marine body is nearly fully enclosed in gas ventilated behind a cavitator. The process of generating a ventilated supercavity, however, is often accompanied by the resonance phenomenon of pulsation. The pulsation of ventilated supercavities, generated by axisymmetric cavitators, was first reported by Song [1] and has been studied extensively by Paryshev [2], among others. Due to his contribution, and his widely cited modeling, pulsating modes of the supercavity are frequently referred to as the Paryshev instability. At extremely high-speeds, cavity gas may be nearly entirely composed of vapor. However, this investigation is restricted to cavities ventilated by a noncondensable gas. This work considers axisymmetric disk cavitators placed in horizontal freestream flows with gravity acting in the vertical direction. Although the interaction of bodies with the cavity closure region is of great interest to researchers and designers (e.g., see Ref. [3]), the cavities studied here do not close on bodies.

For cavities not pulsating, the closure is typically classified as either toroidal vortex shedding (re-entrant jet) or twin vortex. The form of closure determines the main mechanism of gas release from the cavity. Campbell and Hilborne [4] published a frequently referenced paper containing results from a series of supercavitation experiments that support the observation that the product of the Froude number (⁠$Fr$⁠) and the cavitation number (⁠$σ$⁠) can be used, as a rule, to demarcate re-entrant from twin vortex cavities. In this case, the chosen Froude number is based on the freestream velocity and the cavitator diameter [$Fr=U∞/gDn]$⁠, and the ventilated cavitation number is based on the cavity gas pressure $[σ=(p∞−pcav)/((ρ/2)U∞2)]$⁠. Following the standard notation for supercavitation, $U∞$ is the freestream (computational fluid dynamics (CFD) domain inlet) velocity, $Dn$ is the (disk of equivalent drag) cavitator diameter, $p∞$ is the freestream (CFD domain outlet) pressure, $g$ is gravitational acceleration, and $ρ$ is the liquid density. Small values of $σFr$ indicate twin vortex closure and larger ones indicate re-entrant closure. While reducing the ventilation rate in transition from a twin vortex closure to a re-entrant jet closure, the demarcating value was one (⁠$σFr$  = 1). Pulsation is considered a separate regime from the re-entrant or twin vortex closure [5]. When pulsation occurs, it can be during transition from re-entrant to twin vortex or vice versa. The pulsating cavity typically has features of one or both forms of the nonpulsating closures; the distinction being that when cavity pulsation occurs, overall cavity dimensions (i.e., length and diameter) change in a regular, periodic fashion, and waves coinciding with the pulsation period travel axially along the cavity. Pulsation exhibits distinct modal conditions characterized by the number of streamwise pulsation wavelengths appearing on the cavity surface upstream of the cavity closure. Thus, a mode one cavity exhibits a single waveform upstream of the cavity closure, a mode two cavity contains two complete waveforms, etc. Due to the ability of gas to leak steadily from the cavity, pulsation appears to be deterred by the presence of vortex tubes in twin vortex closure; additionally, pulsation may be disturbed by irregular penetration of the re-entrant jet. Thus, pulsation will most likely occur when neither characteristic is too prominent.

Michel [5], in review of the dynamics of supercavities formed behind two-dimensional (2D) wedge and hydrofoil shapes, iterated many of the questions that still are frequently raised regarding ventilated cavity closure modes and cavity pulsation. Even though only reviewing cavity phenomena created by approximately 2D forms, the cavity closures were analogous to those behind axisymmetric cavitators. As was also observed by Campbell and Hilborne [4], Michel discussed the hysteresis regime between the 2D (re-entrant) and three-dimensional (twin vortex) closure modes. Consider the meaning of hysteresis for a supercavity at a fixed Froude condition: As the ventilation rate of the re-entrant cavity is gradually increased, the observed transition to twin vortex closure occurs at a higher flow rate than when starting from a twin vortex condition, and the ventilation rate is similarly decreased. Thus, the mode of cavity closure retains a history and is not merely dependent on ventilation flow rate. Michel stated that when cavity pulsation occurs, it is during transition between 2D and three-dimensional closure.

The linearized modeling of Paryshev [2] results in an eigenvalue problem identifying unstable modes for a given cavity condition. The modeling assumptions required to derive such a modal solution include waveforms traveling downstream from the cavitator detachment location (independent cavity expansions), a compressible gas filled cavity, and some form of gas leakage. The cavity pressure fluctuates with time. It is clear that the model physics should be equally well supported by a multiphase mixture Navier–Stokes model finite-volume approach with a compressible gas component.

Recently, Skidmore et al. [6] conducted an experimental investigation of cavity pulsation. His work focused on the noise generated by various closure regimes for supercavities generated behind a disk cavitator in a closed-loop water tunnel. Skidmore observed that pulsating supercavities exhibit traits of a monopole, with the radiated sound pressure spectrum level generated by pulsation being approximately four orders of magnitude (i.e., 40 dB) greater than comparable twin vortex and re-entrant jet cavities. Additionally, Skidmore provided both experimental conditions and a cavitator geometry that produced the pulsation closure regime. Skidmore [7] observed effects of tunnel size on cavity pulsation and closure mode and also investigated the effects of a partially filled tunnel (i.e., the presence of a free surface) on cavity pulsation. Skidmore observed that, in some agreement with the work of Song [1], the presence of a free surface and tunnel size had a significant effect on pulsation behavior. Song had identified the parameter $β≡[(p∞−pcav)/(p∞−pvap)]$ as a primary indicator of tendency for pulsation. Note that p_vap is the vapor (i.e., partial) pressure and p_cav is the mean cavity pressure. In the more recent work of others [2], this dependence is replaced with that on the nearly equivalent Euler number [$(p∞−pcav)/p∞$].

Most recently, based on experiments conducted in a closed-loop water tunnel facility, Karn et al. [8] presented an investigation of the physical mechanisms responsible for supercavity closure modes. Ventilated supercavity closure modes were classified in eight different regimes. In addition, a foamy cavity regime, i.e., a bubbly mixture (or free jet) not identifiable as a supercavity was considered. Results of pure pulsation closure were not included. However, results at low $Fr$ and highest ventilation flow were classified as pulsating twin vortex. Thus, at the highest supercavity-sustainable ventilation flow (at a ventilation rate higher than pure twin vortex), pulsation was observed. It is interesting to note that, regardless of the ventilation rate, at each $Fr$ condition, all supercavity results fell on a single, corresponding $σ$ value.

The work presented here is intended to demonstrate the capability of the aforementioned numerical physics to capture, in a reasonable way, the different cavity closure modes, and in particular, the previously unstudied pulsation mode. First, the CFD modeling is applied to an idealized, isolated disk cavitator, sufficient to capture the salient physics of the well-studied cavity modes of re-entrant, pulsation, and twin vortex. This is significant because, based on our review of publications, using a finite volume CFD approach, the phenomenon of supercavity pulsation has never, prior to this effort, been reported. (Certainly finite-volume Navier–Stokes CFD is in current and frequent use for industrial high-Reynolds number hydrodynamics.) The modes of cavity closure are demonstrated as are the dependencies of the CFD model to the similarity parameters of $σ$ and $Fr$⁠. The necessary physics are interrogated in the context of the CFD model and the results are evaluated against a linear stability criterion derived by Paryshev. Finally, a specific case studied experimentally by Skidmore is reproduced with near complete similarity.

The investigative method presented here is purely numerical and based on the finite volume method of solving the equations governing multiphase flows. Results are emphasized, and experimental and other modeling data are included for comparison purposes, where appropriate. Computation of incompressible, ventilated supercavity flows using multiphase mixture finite-volume-based fluid mechanics has been developed and published in the literature since the late 1990s and pursued and reported by many including Kunz et al. [9]. A single-fluid, locally homogeneous assumption is made. Therefore, a single pressure and velocity field are shared, at any point in time and space, by all phases in an Amagat-rule mixture model, and no slip is modeled between phases. The Dalton law pressures and densities coexist with Amagat, but are not used or referred to in this work. Similarly, Venkateswaran et al. [10] formulated a compressible model using the density-based mixture model approach. Senocak and Shyy [11] solved a similar set of mixture equations using a pressure-based approach. For single-phase flows, Venkateswaran and Merkle [12] have shown the equivalence of the pressure-based and the preconditioned density-based approaches. The extension to the multiple phase mixture is self-evident. The multiphase mixture model approach to gas liquid flows and large-scale cavitation has proliferated and is found in commercially available and supported software tools such as star-ccm+ [13].

A point of potential interest regards the modeling of turbulence. Kinzel [14] has demonstrated that finite volume-based supercavity modeling is quite sensitive to the turbulence modeling approach and that Reynolds averaging is typically not appropriate for unsteady cavity behavior. In the case of pulsating cavities, at the scales of interest, the physics are unsteady and not appropriately modeled by Reynolds averaging. Wall-bounded turbulent shear layers are not pertinent to the problem of isolated supercavities generated by disk cavitators. Large unsteady features at cavity closure and waves along the cavity, however, are important and must be resolved. As may be shown with mesh and time resolution studies, these are feasibly resolved with nominally available computational resources. Unresolvable, dissipative scales should only be significant downstream from the closure of the main cavity and are, at best, secondary to the cavity dynamics; thus, they are ignored. Although an unsteady turbulence simulation approach, such as large eddy simulation, might be employed, it would require much greater computational resources (e.g., see Ref. [15]) and is not expected to improve the modeling of cavity pulsation. Thus, no proper turbulence model was employed and those resolved scales that may appear turbulent are not rigorously supported with a physical dissipation model. Similar to turbulent dissipation, surface tension quickly becomes significant as the cavity breaks up into generations of bubbly structures. However, Weber numbers characterizing the main cavity are in the range of 1000 or greater. Thus, the cavity interface has a large local radius of curvature and a shape that is determined purely by inertial effects. Surface tension should not be important until well downstream of the point of cavity closure. Of course, this conjecture is fully supported by the presented results.

The commercially available software star-ccm+ was applied to obtain the results presented here [13]. Computational modeling of the gas–liquid homogeneous system is achieved by a second-order accurate space and time finite volume discretization. The transport equations governing fluid motion solved are the mixture continuity equation, the mixture momentum equation, and the liquid mass conservation equation (in volume fraction transport form). Thus, the solution variables are pressure, velocity, and liquid volume fraction. The solution is isothermal, trivializing the energy equation. Thus, for all solutions presented, the temperature throughout the domain is 300 K, and the pure gas sound speed is 347 m/s. The gas is treated as compressible and the liquid as incompressible. Thus, the gas density is determined from a linear relation, depending only on pressure, with a value at 101.3 kPa of 1.18 kg/m³ and the liquid density is constant (997.561 kg/m³), providing the needed closure. As is well known [10], this resulting gas–liquid mixture is extremely compressible and the effective instantaneous sound speed may be derived from the eigenvalues of the mixture equations. Molecular viscosities of water (8.9 × 10⁻⁴Pa·s) and air (1.86 × 10⁻⁵Pa·s) are applied in the governing equations, but for the typical high Reynolds number, supercavitation flow their values have little effect on the results.

As applied here, the gas–liquid interface sharpness is enhanced by modifications of the volume fraction flux term. This approach to gas–liquid interface transport is the volume of fluid (VOF) method. The VOF method, as applied in star-ccm+, consists of the high resolution interface capturing (HRIC) volume fraction convection model and angle corrections [16]. The VOF/HRIC method does not introduce a model of slip velocity or nonequilibrium interface dynamics. If applied properly, the VOF/HRIC method does improve interface sharpness (relative to a standard second- or third-order central or upwind difference advection form). The intent of this work is not to investigate a particular VOF method. However, a few details, regarding the use of the HRIC scheme, are provided to aid in evaluation of the results. The solution maximum convective Courant number was kept below the selected low Courant cutoff value (CFL_l = 0.5), the angle correction value was set to 0.05, and no added compression was used. Otherwise, the guidelines found in the star-ccm+ manual regarding time stepping, and other parameter settings of the VOF scheme were followed [13].

The cavities modeled were ventilated and the minimum absolute pressures found during pulsation integration (assuming nominal atmospheric filled tunnel conditions) were far greater than vapor pressure (CFD cavity minimum pressure was roughly 98 kPa while actual vapor pressure is 3.6 kPa); thus, the formation of vapor was not considered. (Regarding actual tunnel conditions, for a typical case [6], the minimum measured absolute cavity pressure was roughly 90 kPa.) The time integration of the segregated, pressure-based system was completed by the standard implicit (subiterative) method [13]. For all results presented here, five subiterations were completed per time step, and, in addition to the care taken to maintain the Courant criterion, the time step convergence was monitored to ensure that at least one order of magnitude reduction in governing equation residuals was achieved per time step. Discretization convergence in space and time was also tested.

The formulation of the volume fraction equation is made with application of the chain rule to transported quantities. Although a mathematically correct operation resulting in a more computationally robust formulation [13], the result may not be algebraically consistent with jump conditions and is therefore called nonconservative. A nonconservative form may be of concern when one or more of the mixture elements are compressible or at a resolved discontinuity in volume fraction. The term of concern arises in the governing mixture mass conservation equation. A source term is present representing quantities that are nonzero when the velocity field is not divergence free [13]. For a nonconservative method, it is expected that some compressible flows, such as those involving strong shocks, may be pathological. Thus, sensitivity to the nonconservative form should be resolved by conducting a mesh sensitivity study. The flows of concern here are of very low freestream (single-phase) Mach number, and shock waves are not expected. In fact, it has been pointed out and will be verified that for ventilated supercavity pulsation, the pressure distribution of gas within the supercavity should be mostly time-varying but, at any one time, nearly spatially uniform. Thus, the physics of interest do not involve strong compression or expansion waves.

Numerical results are presented for both an isolated, canonical disk cavitator and a model based closely on the geometric and fluid similar configuration of an actual water tunnel experiment. As it has been well established that a disk cavitator is a representative and useful benchmark for ventilated supercavitating configurations of interest, the phenomenon of pulsation and the ability of the CFD approach to capture the expected pulsation behavior will first be presented. By adjusting the controlling similarity parameters of supercavitation, and reviewing the results, the phenomenological goodness of the approach is verified. Finally, as a further, more precise validation of the CFD approach, a pulsation result is presented matching most of the similarity parameters to a water tunnel experiment.

For all results presented, similar boundary conditions were applied. The cavitator was set in a straight right circular cylinder tunnel section with inviscid walls (i.e., symmetry boundary condition). The tunnel inlet boundary was a flat, circular, planar section, the physical condition was constant total pressure and a volume fraction of pure liquid, while the numerical velocity at the boundary was directed axially (i.e., normal to the boundary surface). The tunnel outlet was also a flat circular plane section, and the physical condition was constant static pressure. The wetted and lateral facing surfaces of the cavitator were treated as no-slip walls. The aft face of the cavitator was the ventilation gas inflow boundary. The ventilation volume fraction was set to pure gas, and the ventilation flow entered normal to the boundary and the velocity was therefore chosen to match the desired C_Q. Thus, the ventilation volumetric flow rate was constant and determined by the specified inlet velocity boundary condition.

A representative result of each kind of cavity closure, all captured on the same CFD mesh, is given in Fig. 1. Snapshots from re-entrant, twin vortex, and pulsating closure cavity modes are shown from both water tunnel testing and the current computational modeling. The differing CFD results were obtained by varying the ventilation flow rate but maintaining Froude similarity. However, the CFD model and water tunnel at similar closure conditions are not at complete fluid similarity, i.e., do not represent the same Froude, ventilation rate, and cavitation number. This figure conveys the phenomenological capability of the CFD.

The effect of Froude number on pulsation (in the model simulations) was strong. The CFD cases shown in Fig. 1 are all at a Froude number of 12.3. The computational mesh applied in these cases was relatively coarse. However, it is clear that the flow details characteristic of each type of cavity are captured. It is expected that the rate of gas leakage and the fine structure of the cavity closure may be better resolved with a finer mesh. The results in this figure seem to indicate that the current CFD-based modeling is physically sufficient to resolve all three modes of cavity closure.

Computational results were obtained for supercavities formed behind a circular disk cavitator placed in a flowing fluid domain. No turbulent boundary layer-resolving mesh was included on the cavitator itself, as the boundary layer flow on the front side of a disk cavitator is always in a favorable pressure gradient and should be laminar. A mesh sensitivity study was conducted using the standard Richardson extrapolation-based approach for a representative Froude 24.5 case. Due to the similar physics in all presented cases, findings of this Froude 24.5 mesh study are intended to support the numerical convergence of all other results.

In Fig. 2, a series of three computational meshes used during the investigation are shown. The mesh is shown on the cavitator, inflow and outflow boundaries, and to best illustrate the mesh density in the region of interest, the mesh intersection with a plane bisecting the domain is shown. No physical symmetry assumptions were made; the bisecting plane is used merely to present the mesh density. Both the entire extent of the domain and a close-up of the area requiring extra mesh resolution are shown. In this case, the computational tunnel is roughly 53 times larger than the cavitator diameter.

The three computational meshes shown in Fig. 2 are sequenced with a factor of two refinement appropriate for application of a Richardson extrapolation to estimate the numerical error. The grid convergence index approach credited to Oberkampf and Roy (GCI-OR) and found in Phillips and Roy [17] is applied to estimate the fine mesh error and uncertainty. For the series of solutions, the inlet boundary condition was held at a constant stagnation pressure, relative to the outlet static pressure, to yield a mean inlet velocity of 15 m/s. The cavitator diameter was 3.81 cm, resulting in a Froude number of 24.5. The mean cavity pressure, as measured at a numerical probe located 0.67 diameters downstream of the cavitator face, was used as the convergence measure. In Table 1, the nominal mesh size Δh is reported along with the time step size Δt, the mesh refinement factor r_ij, the observed cavity closure, and the mean cavity pressure (in terms of σ). For the mesh strategy used, the reported Δh was a controlling global dimension, and cell dimensions throughout the mesh were scaled on that parameter. Using a direct application of Richardson's extrapolation method, the observed order of convergence, $p̂=6.31$⁠. The large value of convergence exponent is due to the coarse mesh being outside of a reasonable asymptotic range. The numerical method is formally second-order accurate. Therefore, following Ref. [17], the error estimation exponent is upper limited to 2 and the inflated factor of safety, equal to 3, is applied, indicating a fine mesh error of 0.2% with uncertainty of 0.6% (both of dynamic pressure).

No matter the attempted initial condition or time step size, the coarse grid integration was insufficient to capture pulsating supercavities, resulting in a model re-entrant condition. The difference between the fine and coarse mesh histories is clear and shown in Fig. 3. The coarse mesh pressure history is at a far lower mean level and has a very different (not periodic) unsteady content from the medium or fine mesh results. The medium mesh results exhibit pulsation at 19.94 Hz, and the fine mesh results pulse at 20.14 Hz. Finally, note that time step size was studied for convergence on the fine mesh and also was sufficient based on requirements of the as-applied HRIC scheme (maximum convective Courant number less than 0.5).

Traditional modeling of supercavities employs the concept of independent expansions. Typically, this approach begins with the assumption that at a given moment in time, cavity pressure is nearly uniform throughout the cavity [18]. For a large, continuous gas-filled region with low Mach number everywhere, this is intuitive. This notion is by no means a requirement for a CFD approach based on finite volume or other similar discretized transport equation-based approach to capturing pulsating supercavities. However, if strong compression and expansion wave forms are present, such as one might actually find in a frothy, dynamic, liquid–gas mixture, the nonconservative formulation applied in obtaining the results presented here might merit further scrutiny. To illustrate the nature of the pressure field in a gas-filled supercavity, consider Fig. 4. Here, the pressure histories at six probe locations distributed along the length of the $Fr=24.5$ pulsating cavity are presented. These records clearly indicate that the cavity pressure is oscillating as a roughly lumped parameter, dependent only on time, and that discrete or strong nonlinear pressure waveforms do not play a role in cavity pulsation. Thus, as anticipated, within the supercavity, at any moment in time, there are negligible spatial pressure gradients. Although traces at five of the six probes are nearly indistinguishable, there are slight differences observable in the pressure history of probe located in the furthest aft part of the cavity, i.e., in the two-phase closure region. These computational observations are consistent with the traditional supercavity modeling theories such as those developed by Logvinovich [18] and Paryshev [2].

Sample time histories of simulated cavity pressure during twin vortex, pulsation, and re-entrant cavity closure are given in Fig. 5. The cavity pressure was determined with a numerical probe placed 0.67 cavitator diameters downstream of the front face of the cavitator. Note that the numerical cavitator itself was a finite thickness disk. Based on the results in Fig. 4, only a single probe aft of the cavitator was deemed necessary. The behavior of the simulated pressures shown in Fig. 5 appears to be similar in nature to those measured experimentally by Skidmore et al. [6], also shown in the figure. The twin vortex cavities for both simulations and experiments have a roughly constant pressure coefficient. The re-entrant cavities for both simulations and experiments have some spectral content and noise, but are neither regular nor sinusoidal. The pulsating cavities for both simulations and experiments are distinctly regular and periodic.

An overview of closure modes for a range of cases is given in Fig. 6. The same results are also presented in Table 2. Computationally, at least, Froude number influences the observed results. At the lowest Froude number, 5.78, no hysteresis or pulsation was observed, and cavity transition from re-entrant to twin vortex coincides with, approximately, $σFr=1$⁠. At the middle Froude number, 12.3, some hysteresis was observed and the transition between re-entrant to twin vortex closure also coincides with $σFr=1$⁠. At the highest Froude number, 24.5, pronounced hysteresis and pulsation were observed. The transition to twin vortex occurs close to $σFr=1.2$ while the return to re-entrant is closer to $σFr=1$⁠. The Froude 24.5 hysteresis loop is highlighted in Fig. 6. Thus, the highest Froude number case agrees most completely with the findings of Campbell and Hilborne [4]. To determine the demarcation (⁠$σFr$⁠) value, pulsation cases were neglected. Thus, as expected, the pulsation cases fall in the hysteresis region. Many of the pulsation cases appeared to also be re-entrant in nature, while a few appeared to have buoyant, twin vortical structures at closure. It is believed that the difficulty in capturing pulsation at the low Froude condition is a mesh resolution issue. However, this was not resolved.

An image from a CFD simulation representing a typical twin vortex cavity is shown in Fig. 7. The mesh density on the cavity surface is illustrated in the figure. The ridge formation along the axial axis of the supercavity on the lower surface can be clearly observed. This ridge can be observed to grow in height until reaching the upper surface at the closure of the supercavity, which is both expected and observed experimentally. Termination of the captured vortical structures in the figure is due to reduced mesh density further downstream. Hence, for the cavity lengths considered, the mesh density was adequate to capture all supercavity physics of concern.

A review of the modeling of Paryshev [2] suggests that cavity compressibility is fundamental to pulsation. However, it is interesting to investigate the possibility that pulsation is not dependent on the cavity gas acting as a compressible spring. The CFD model is used to test this conjecture. In Fig. 8, snapshots from a case that was initially pulsing with compressible physics are presented. When restarted with incompressible physics, the modeled cavity quickly transitioned to re-entrant behavior with no pulsation. A snapshot with incompressible physics and re-entrant behavior is shown next to a snapshot of the solution with compressible physics. The incompressible cavity surface waves are irregular and nonperiodic. The compressible gas case exhibits clear second-mode waves. Also, in the figure, cavity pressure histories from the incompressible and compressible integrations are shown. Clearly, the compressible history appears regular and periodic (it is pulsating), and the incompressible history is not regular (it is twin vortex). Although not an exhaustive proof that pulsation is a compressible phenomenon, it is compelling.

A series of snapshots from a CFD simulation representing a second mode pulsation condition is shown in Fig. 9. This cavity, at any point during the simulation, contains two complete traveling waves along its length, upstream of closure (exhibits mode two or second-order pulsation). In this case, a stability criterion presented by Paryshev [2] was also evaluated. The Paryshev linear stability criterion at the conditions shown in Fig. 9 results in a scaled ventilation coefficient: $(CQσ/CD)κ=0.160$⁠. Here, C_Q is the dimensionless ventilation coefficient, C_D is the cavitator drag coefficient, and σ is, again, the cavitation number. Values less than 0.0160 are expected to pulsate. The CFD value of $CQσ/CD$ at the condition shown is 0.0157, as indicated in Fig. 10. Considering the approximate nature of the Paryshev derivation, this value is essentially on the stability boundary. Interestingly, at a ventilation rate just 4% lower, also shown in Fig. 10, the CFD solution transitions to re-entrant jet closure. This re-entrant jet cavity has a value of $CQσ/CD$  = 0.0255, well above the Paryshev threshold for stability. Thus, the CFD results are in concurrence with the Paryshev criterion.

A primary computational objective was to capture supercavity pulsation at water tunnel fluid similarity conditions. The physical model was based on a previously published experiment [6]. In Fig. 11, a schematic of the test apparatus is presented. The apparatus was installed in a 0.305 m water tunnel. To capture this computationally, a 3.9 million cell mesh, shown in Fig. 12, was applied. The CFD model included geometrically similar representation of the sting and cavitator from the experiments, but neglected the strut. In addition, the length of the simulated tunnel test section was greater, neglecting the contraction and expansion regions, and the diameter was also greater. The computational tunnel diameter was 0.61 m. Tunnel walls were treated as inviscid. Thus, relative to the physical experiment, tunnel blockage was neglected. Finally, the details of the actual ventilation ports were neglected and gas flow was introduced into the CFD domain at a uniform volumetric rate at the cavity pressure from the back planar face of the cavitator as indicated in Fig. 12 (based on the boundary condition formulation [13]).

A photograph from the experiment, with a pulsating cavity, is shown in Fig. 13. The cavity was generated at a Froude number of 4.5, both following the methodology and based on the instrumentation of Ref. [6]. The cavity is mode two based on the observed presence of two waves on the cavity interface upstream of the cavity terminus. The terminus is defined as the location where the apparent cavity diameter becomes less than that of the cavitator. The cavitation number was measured to be σ = 0.31. This measured value is in rough agreement with the correlation based on cavity diameter ([19], Eqs. (3)–(5)), σ_May = 0.36. Matching the Froude number to experimental condition, a computational pulsating supercavity was generated and it is also shown in Fig. 13. The CFD cavity interface is visualized using both isosurfaces of gas volume fraction equal to 0.5 and to 0.25. The 0.25 isosurface is presented to more clearly illustrate the mode 2 cavity shape (similar to the experiment). This CFD result is at σ = 0.26. Thus, the CFD mean pressure is in rough agreement with the mean pressure determined in the experiment.

In Fig. 14, the experimental and computational time histories of internal cavity pressure are presented. Clearly, both histories are dominated by a single common frequency, and the amplitudes of oscillation are approximately the same. In Fig. 15, and all reporting here, the pressure spectra are presented in dB referenced to 1 μPa. The internal cavity pressure as well as the near-field radiated noise spectra is shown. All signals indicate pulsation at a frequency of 38 Hz. In order to ensure that both the experimental and computational signal underwent the same spectral analysis, an identical length of time record was taken of all pressure signals. These were then up-sampled to match the time step size of the simulation, 1 μs. This results in the cavity interior pressure spectrum level were 176.1 dB at the pulsation frequency, while the pressure spectrum level measured by the window hydrophone, also at the pulsation frequency, was 156.5 dB. The 19.6 dB reduction in spectrum level at the window hydrophone location is due to spherical spreading of the sound waves from the cavity interface [6].

Much like the experimental cavity, the CFD cavity pressure, also in Fig. 15, measured via a pressure probe located 0.63 diameters aft of the cavitator forward face, was dominated by a single frequency oscillation. The corresponding spectrum of the cavity interior pressure computed from the CFD simulation time series indicates pulsation at a frequency of 40 ± 1 Hz, roughly similar to the experimental 38 ± 1 Hz measurement. The spectrum level of the simulated cavity interior pressure at the pulsation frequency was found to be 177.3 dB, close to the measured spectrum level of 176.1 dB. Two approaches were taken to estimate the pressure spectrum level at the window hydrophone location; one accounting for spherical spreading of the sound waves from the cavity interface and the other evaluating Eq. (5) of Ref. [6] for a monopole source with the simulated cavity volume velocity. For the former method, the mean radius of the simulated cavity is 1.69 cm and the window hydrophone is 22.4 cm away from the cavity, giving an estimated pressure spectrum level of 154.9 dB at the hydrophone location. This is in reasonable agreement with the measured 156.5 dB (Fig. 15). Applying the method of Ref. [6] to the CFD simulation, the cavity radius varies by ±0.61 cm about the mean value of 1.69 cm at 40 Hz yielding a volume velocity of 0.0019 m³/s and a pressure spectrum level of 158.6 dB (Fig. 15). The volume velocity prediction agrees well with both the level predicted from the simulated cavity interior pressure and with the measured level.

However, there is one similarity parameter that is not properly matched between computation and experiment. The computational ventilation rate is much lower than that required in the experiment, C_Q-CFD = 0.09 and C_Q-EXP = 0.91. Thus, there are small discrepancies in unsteady cavity pressure (roughly 5% of mean dynamic pressure with negligible differences in frequency), and apparent cavity size and shape, but extreme differences between computational and water tunnel ventilation requirements. In the water tunnel testing, the tunnel was at near atmospheric pressure and the ventilation gas, and air was injected at near room pressure and temperature. The CFD model gas is injected at the cavity pressure (as is required by the velocity inlet boundary condition), and obeys the linear, isothermal state relation described in the “Computational Modeling” section. The actual water tunnel injection orifices are illustrated in Fig. 11, and the CFD model injection is shown as the entire back face of the cavitator (Fig. 12). Therefore, for a given volumetric flow rate, the CFD injected gas must be introduced to the cavity at a much lower velocity than in the actual testing.

The supercavity ventilation requirement (for a given cavitation number) is a physical quantity that has proven difficult to capture with finite volume CFD. It is of interest to understand the mechanisms affecting the gas entrainment rate as well as to determine reliable correlations to predict it. Investigations of cavity gas entrainment mechanisms have been published by Spurk [20] and Kinzel [21]. It is suspected that the air entrainment rate from a ventilated cavity is dependent on many physical conditions beyond the type of cavity closure. Some are difficult to computationally capture but important features are the presence of small disturbances and waves on the cavity air–water interface (i.e., cavity roughness). In addition, although it may be computationally feasible to capture the largest structures of cavity breakup, many smaller bubbles frequently form in the closure region and elsewhere around a real cavity accounting for added gas loss. Finally, it must be noted that the air entrainment rate required for a given cavitation and Froude number will also be somewhat dependent on more difficult to control physical features such as water quality and freestream turbulence.

The three main dynamic closure modes (i.e., re-entrant jet, pulsation, and twin vortex) of supercavities developed behind an axisymmetric cavitator may be properly resolved using finite volume CFD. The physics required to capture the re-entrant jet and twin vortex conditions are incompressible and have been demonstrated previously [22]. The physics required to capture pulsation are compressible (in the gas phase). This is in agreement with the basic theory of Ref. [2] and has been demonstrated here. A number of unsteady supercavity computational simulations have been completed and mesh convergence has been demonstrated. The conditions were determined by varying the Froude number (based on upstream velocity and cavitator diameter) and the ventilation flow rate. The CFD solutions presented here collectively fit the expected interdependencies of cavitation number, ventilation rate, and Froude number (Fig. 6). Three different Froude conditions were modeled, and for conditions exhibiting pulsation, re-entrant jet-twin vortex hysteresis was observed.

The computed dynamics of pulsation have been correlated between cavity shape and pressure time histories. As confirmation supporting the applicability and accuracy of the chosen CFD approach, the conditions of cavity closure, including unsteady pulsation, were successfully and pleasingly captured using cavity size and shape histories and pressure histories. In addition, by assemblage of unsteady results in C_Q-σ phase space, the expected patterns [4] were observed. Similarly, the results were found to be in agreement to a linear stability criterion due to Paryshev [2].

Finally, a computationally simulated ventilated, pulsing cavity was compared to an experimentally similar pulsing cavity at identical Froude conditions. Computational ventilation rate was adjusted to capture cavity pulsation, and the results were compared. The cavitation number, cavity pulsation order, pulsation frequency, cavity interior pressure, and radiated sound pressure spectrum levels were all found to be captured well as represented in the computational results. On the other hand, the computational air entrainment rate (ventilation rate) differed greatly between CFD and the experiment. It is conjectured here that this discrepancy is due to the many real physical mechanisms of cavity gas entrainment that are extremely difficult (if possible) to simultaneously resolve in a finite volume CFD method. Fortunately, in order to properly capture all of the modes of cavity closure, the expected dependencies of these modes on the ventilation rate, and also capture the pulsation phenomenon, it seems only necessary to capture some of these air entrainment mechanisms. Thus, there are mechanisms of gas leakage from the supercavity that have only a tertiary effect on cavity closure mode and dynamics, and so long as a controlling mode of gas leakage is sufficiently captured, the cavity dynamics may be resolved in a finite volume approach.

It has been shown that this method of computationally simulating pulsating supercavities is capable of modeling the majority of dynamics and closure of experimental supercavities; thus, considering the limits of what may be done in water tunnels and tow basins, the present CFD method may be applied to investigate physical supercavity dynamics beyond what is experimentally feasible.

Future work should attempt to improve the understanding of the sensitivity of closure mode and cavity dynamics to water tunnel blockage and the presence of a free surface. Additional efforts should be made to better model air entrainment mechanisms in a finite volume approach, perhaps via an eddy viscosity or similar concept. Finally, it should be straightforward, using the current approach, to investigate control of pulsation.

This work was sponsored by the Office of Naval Research under Grant No. 00014-11-10322 and monitored by Dr. Ronald D. Joslin.

$CD$ =

coefficient of drag [$Fx/((ρ/2)U∞2π(Dn2/4))]$

$Cp$ =

pressure coefficient [$(p−p∞)/((ρ/2)U∞2)]$

$CQ$ =

dimensionless ventilation rate [$Q/(U∞Dn2)]$

$Dn$ =

cavitator diameter (diameter of a disk of equivalent drag)

$Fx$ =

cavitator drag

$Fr$ =

Froude number [$U∞/gDn]$

$g$ =

gravitational acceleration

$pcav$ =

cavity pressure

$p∞$ =

freestream pressure

$Q$ =

volume flow rate of ventilation gas

$U∞$ =

freestream velocity

$ρ$ =

water mass density

$σ$ =

cavitation number [$(p∞−pcav)/((ρ/2)U∞2)]$