## Abstract

We conducted a numerical simulation of ventilated supercavitation from a forward-facing cavitator in unsteady flows generated by a gust generator under different gust angles of attack and gust frequencies. The numerical method is validated through the experimental results under specific steady and unsteady conditions. It is shown that the simulation can capture the degree of cavity shape fluctuation and internal pressure variation in a gust cycle. Specifically, the cavity centerline shows periodic wavelike undulation with a maximum amplitude matching that of the incoming flow perturbation. The cavity internal pressure also fluctuates periodically, causing the corresponding change of difference between internal and external pressure across the closure that leads to the closure mode change in a gust cycle. In addition, the simulation captures the variation of cavity internal flow, particularly the development internal flow boundary layer along the cavitator mounting strut, upon the incoming flow perturbation, correlating with cavity deformation and closure mode variation. With increasing angle of attack, the cavity exhibits augmented wavelike undulation and pressure fluctuation. As the wavelength of the flow perturbation approaches the cavity length with increasing gust frequency, the cavity experiences stronger wavelike undulation and internal pressure fluctuation but reduced cavitation number variation.

## 1 Introduction

Supercavitation is a special case of cavitation in which a cavity is large enough to encompass the object traveling in the liquid [1,2]. It provides a promising approach to augment speed of underwater vehicles with drag reduced by as much as 90% [3]. A supercavity can be formed at low speeds by ventilating noncondensable gas into the low-pressure region around a cavitating object (i.e., cavitator), termed as artificial or ventilated supercavitation [4]. Ventilated supercavitation is generally characterized by using several nondimensional parameters, including ventilated cavitation number $\sigma =2(p\u221e\u2009\u2212p)/(\rho wU\u221e2)$, Froude number $Fr=U\u221e/gdc$ and air entrainment coefficient $CQ=Q\u02d9/U\u221edc2$, where $p\u221e\u2009$ is the ambient pressure upstream of the cavitator, $pc$ is the cavity pressure measured inside the cavity, and $\rho w$, $U\u221e$, and $g$ refer to the water density, the upstream incoming flow velocity at the test section, and gravitational acceleration, respectively, *d*_{c} corresponds to the cavitator diameter, and $Q$ denotes the volumetric ventilation rate. A number of researchers have studied the ventilated supercavitation to derive a relationship between the closure modes and the ventilation demand [5–7]. However, there is no consensus on the analytical models for the supercavity closure and the gas ventilation [8], and no consensus on the formation conditions for different closure modes [4,9]. Karn et al. [10] investigated systematically the closure modes of ventilated supercavity under different flow and ventilation conditions using a backward-facing cavitator model at Saint Anthony Falls Laboratory (SAFL). Their study not only revealed several new closure modes, but also provided a unified physical framework to explain the variation and transition of different closure modes. Using particle image velocimetry, a recent investigation [11] studied the internal flow of a ventilated supercavity under different closure types. The internal flow results suggest that at the upstream of the location of the maximum cavity diameter, the gas enters the forward flow (including the internal boundary layer and the forward moving portion of the ventilation influence region) from the reverse flow, while at the downstream of that location, the gas is stripped from the internal boundary layer and enters the reverse flow due to the increasing adverse pressure gradient in the streamwise direction. Combined with visualization results of cavity geometry and closure patterns, the results given previously are able to explain the influence of the cavity gas leakage mechanisms on cavity growth and closure transition. However, due to immerse technical difficulties, the experimental study of cavity internal flow has yet calculated the pressure variation inside the cavity, and the cavity internal flow under unsteady flow conditions has not been conducted.

In practical applications, the supercavitating vehicle is expected to be operated under unsteady external flow conditions. Particularly, when the underwater vehicle travels near the sea surface, the supercavity may encounter an unsteady incoming flow induced by the sea waves. However, there is only a very limited number of experimental works on supercavitation in unsteady flows. For example, Lee et al. [12] experimentally studied the ventilated supercavitation around a vehicle pitching up and down with the emphasis on the interaction between the vehicle aft body and supercavity boundary. Lee et al. [13] duplicated several sea states using a periodic gust generator in the high-speed cavitation tunnel at SAFL to investigate the effects of unsteady flows on the shape deformation of axisymmetric supercavities. Using the same setup, Karn et al. [14] provided further insight into the dependence of supercavity closure on the flow unsteadiness and found that the incoming unsteady flows would cause the transition of the supercavity closure and its shape. Systematic experiments were performed by Shao et al. [15] to explore different states of a ventilated supercavity under various unsteady flow settings. Although synchronized measurements of cavity shape and pressure are conducted in the past [10,13–15], the pressure inside the ventilated supercavity is measured at a fixed location and the information regarding to the internal velocity and pressure field is still quite scarce in periodic gust flows due to the limitation of diagnostic approaches.

With the advancement of computational fluid dynamics, numerical simulation has been implemented to provide a more detailed understanding of the cavity internal flows and the characteristics of supercavitation under the conditions which cannot be readily achieved by the experiments. Specifically, using Reynolds-averaged Navier–Stokes (RANS) simulation, Cao et al. [16] investigated the pressure distribution inside a ventilated supercavity under different ventilation rates and the blockages in a closed-wall water tunnel. The work highlighted the strong pressure variation near the closure region of the ventilated supercavity, which was potentially connected with different closure modes as pointed out in Karn et al. [14]. Rashidi et al. [17] used the volume of fluid method coupling with Youngs' algorithm in the advection of the free-surface to investigate physics of ventilated cavitation phenomena (i.e., the cavity shape, the gas leakage, and the reentrant jet), and the reentrant jet was found to cause the variations of the cavity pressure and the transient flow behaviors such as the cavity detachment and the internal flow inside the cavity. Using a free surface model and a filter-based approach, Wang et al. [18] numerically investigated the gas leakage behavior and reentrant jet dynamics for a vehicle body from partial cavitation to supercavitation. They found two mechanisms of the gas leakage at different Froude numbers, termed toroidal vortices mode and two hollow tube vortices mode, and the reentrant jet at the cavity tail was shown due to the recirculation of water into the cavity. By using the shear stress transfer *k*–*ω* turbulence model with a mass transfer modeling, the natural cavitation and supercavitation were studied around submarine hull shapes with/without sail and appendages [19]. The sail and appendages were found to promote the supercavity occurrence. Park and Rhee [20] used the full cavitation model with the van Leer scheme to simulate the natural supercavity around a two-dimensional wedge-shaped cavitator for various wedge angle and cavitation numbers, and the results were in good agreement with those obtained with an analytic solution and the potential flow solver. However, the effects of the incoming flow unsteadiness on the internal flow and pressure fields inside the supercavity have not yet been numerically explored and investigated in detail.

This paper presents some numerical results for the ventilated supercavity in a periodic gust flow with different fluctuating amplitudes and frequencies, particularly focusing on the dependence of the supercavity deformation and internal field variations upon the periodic gust flow. The paper is structured as follows: Sec. 2 provides a detailed description of the numerical approach including the governing equations and the simulation setup as well as the experiment for validating the simulation. Subsequently, Sec. 3 presents the results of the supercavity deformation and transient internal flows under various unsteady flow settings. Finally, a summary and discussion of the results are provided in Sec. 4.

## 2 Research Approach

### 2.1 Governing Equations and Computational Setup.

where $p$ is the pressure, $\mu k$ is the molecular dynamic viscosity of phase $k,$ and $\mu t,\u2009k$ corresponds to turbulent eddy viscosity of phase *k* which is given by the standard $k\u2212\epsilon $ turbulence model [21]. $Sk,buoy=(\rho a\u2212\rho w)g$ corresponds to the momentum source due to buoyancy forces. $Mk=CD\rho mAm|ua\u2212uw|(ua\u2212uw)$ is the total interfacial drag forces acting on the phase water due to the presence of phase air, where the nondimensional drag coefficient $CD$ is treated with the default value in ansyscfx of 0.44, $\rho m=\u2211\alpha k\rho k$ is the density of the mixture and the interfacial area per unit volume $Am$ for the free surface model is calculated by $Am=|\u2207\alpha k|$.

The time-dependent supercavitating turbulent flows are simulated using RANS equations together with the standard $k\u2212\epsilon $ turbulence model. The high-resolution scheme [22] is used for the advection terms. This scheme is a second-order scheme, which can locally switch to first-order scheme to prevent numerical oscillations near critical high-density gradient areas [23]. So this scheme is both accurate and robust since it only reduces to the first order near discontinuities. The central difference scheme is used for the diffusion terms and the second-order backward Euler scheme is used for the transient terms. The discrete continuity and momentum equations for the multiphase flow field are solved together without iterations and corrections which improve the stability of the numerical scheme.

As shown in Fig. 1, a forward facing cavitator with the diameter of 20 mm, described by Ref. [4], matching the model used in the experiments is employed to generate the cavity and the three-dimensional computational domain is 1.55 m in length with the same square cross section (0.19 m × 0.19 m) of the experiments. The streamwise and vertical direction corresponds to the *z*-axis and *y*-axis of the Cartesian coordinate system, respectively. The ventilation ports at the rear part of the cavitator are six circumferentially placed holes in the experiment. Such configuration is simplified to an annular band in simulation, which is an effective method to improve the local mesh quality and has been confirmed of its feasibility in supercavitation simulation in Ref. [16]. For the boundary conditions, the velocity components are used at the inlet that the streamwise velocity is $U\u221e=8.5\u2009m/s$ with zero spanwise velocity (i.e., *x*-direction). The unsteadiness of the flow in the current paper is introduced through adding a vertical velocity term $vg(t)$ at the inlet as a function of the pitching angle of the gust generator fitted to the experimental data [13,24]. No-slip conditions are employed on the other four boundaries of the computational domain except the inlet and outlet. A detailed description of the gust generator is provided in Sec. 2.2. The pressure is adjusted at the outlet to match the cavitation number in experiments. The air flowrate ($\u2009CQ=0.15$) is set at the ventilation port of the cavitator for all the cases. Since the vertical velocity term is updated every time-step during the transient simulation, it is necessary to achieve the convergence in each time-step with a criterion that the calculation for each time-step is 1000 iterations with the root-mean-square (RMS) residuals below 10^{−5}.

Structured grids are generated in the computational domain as shown in Fig. 2(a). The scalable wall function is used herein with a requirement of 11.06 < *y*^{+} <300, where *y*^{+} is the dimensionless distance from the solid wall of the cavitator, so grids around the cavitator are refined to allow a sufficient resolution to resolve the detailed physics of the supercavity.

In order to determine an optimal grid resolution, three cases are tested under the flow conditions of $Fr=18.7$, $\u2009CQ=0.15$, $dc=20\u2009mm$, and $\sigma steady=0.20$. Three grid quantities are determined with a constant grid refinement ratio $r=1.3$ in all three directions according to Table 1. As shown in Fig. 2(b), cavity maximum diameter locations (point 1 and 2) and a point near the cavity closure region (point 3) are chosen as the monitoring points for the comparison of results from different cases. The grid convergence index (GCI) of velocity and static pressure at monitoring points [25] is introduced to estimate the uncertainty. As shown in Table 1, the uncertainty estimated by GCI method along the streamwise direction with has a value less than 5% which demonstrates that the simulation results are almost independent of the grid resolution. Additionally, the comparison results of cavity shape (defined the contours of air volume fraction equals to 0.1) under three grid quantities are closely matched with each other with small differences at the cavity closure region (Fig. 2(b)). Note that the further refinement of grids will lead to substantial increasing of computational cost and the instability of the numerical simulation. Therefore, we choose the case 2 as the final mesh in the simulation corresponding to about 1.6 × 10^{6} elements.

Mesh | Grid quantities | u_{1} | p_{1} | u_{2} | p_{2} | u_{3} | p_{3} |
---|---|---|---|---|---|---|---|

Case 1 | 2,644,817 | 9.04 | 58203.00 | 9.33 | 58305.60 | 8.61 | 58245.90 |

Case 2 | 1,555,775 | 9.03 | 58484.60 | 9.25 | 58613.00 | 8.74 | 58537.30 |

Case 3 | 915,161 | 9.01 | 58533.30 | 9.25 | 58674.10 | 8.66 | 58669.60 |

GCI | 1.29% | 0.13% | 0.02% | 0.16% | 3.67% | 0.52% |

Mesh | Grid quantities | u_{1} | p_{1} | u_{2} | p_{2} | u_{3} | p_{3} |
---|---|---|---|---|---|---|---|

Case 1 | 2,644,817 | 9.04 | 58203.00 | 9.33 | 58305.60 | 8.61 | 58245.90 |

Case 2 | 1,555,775 | 9.03 | 58484.60 | 9.25 | 58613.00 | 8.74 | 58537.30 |

Case 3 | 915,161 | 9.01 | 58533.30 | 9.25 | 58674.10 | 8.66 | 58669.60 |

GCI | 1.29% | 0.13% | 0.02% | 0.16% | 3.67% | 0.52% |

Note: Point 1 and point 2 are at the maximum-diameter location and point 3 is near the cavity closure. The units for velocity and pressure are m/s and Pa, respectively.

### 2.2 Experimental Setup.

The validation experiments are conducted in a high-speed water tunnel in the Saint Anthony Falls Laboratory at University of Minnesota with details provided in the previous studies [10,14]. It is a closed recirculating facility and the length of the horizontal test section is 1.20 m long and the same cross section with the simulation. A schematic of the experimental setup is shown in Fig. 3. Two oscillating NACA0020 hydrofoils with chord length of 40 mm are placed 67 mm upstream of the cavitator to generate a periodical varying unsteady flow. The hydrofoils can oscillate in phase driven by a pivot arm which is connected to an eccentricity flywheel outside the tunnel. It is worth noting that unlike the simulation, the flapping motion of the hydrofoils provides a test section pressure variation associated with varying blockage ratios of the hydrofoils besides the vertical velocity term. Instantaneous velocities of the periodic gust flows are measured by laser Doppler velocimetry in Refs. [13] and [24] with the vertical velocity component fitted to a sinusoidal function, i.e., $vg(t)=vgmaxsin(2\pi fgt)$, where $vgmax$ is the maximum vertical velocity (i.e., amplitude) measured near the leading edge of the stationary hydrofoils controlled by the angle of attack ($\gamma g$) of two hydrofoils, and $fg$ is the frequency of periodic flows (i.e., wavelength) corresponding to the rotational speed of the motor. A forward facing model that matches that of the simulation is employed in experiments to generate a ventilated supercavity with smooth surface. The air flow is injected from the ventilation ports at the cavitator and controlled by a mass flow controller. During the experiments, two sensors are employed to measure the pressure before the gust generator and within the cavity 64 mm behind the cavitator. The pressure sensors are sampled at 1000 Hz with the uncertainty around ±100 Pa. The high-speed imaging is synchronized with the pressure measurements through the LabVIEW software. Detailed experimental setup and data acquisition procedures can be found in Ref. [15].

### 2.3 Validation of the Numerical Approach.

First, the ventilated supercavity in the steady incoming flow is used to validate the present numerical approach. The conditions are $Fr=18.7$, $\u2009CQ=0.15$, $dc=20\u2009mm$ and $\sigma steady=0.20$ measured in experiments, under which the cavity has a twin-vortex closure type. The numerical simulation is carried out with the same boundary conditions as the experiment by adjusting the outlet pressure to match the cavitation number. The experimental and simulated supercavity geometries are compared in Fig. 4 with two geometrical parameters proposed by Brennen [26], i.e., maximum diameter ($Dmax$) and half-length (*L*_{1/2}) of the supercavity which is a horizontal distance between the cavitator and $Dmax$. It is shown that the simulated cavity has maximum diameter of 61 mm and half-length of 200 mm, while the experimental results are 62 mm and 204 mm, respectively. The predicted supercavity shifts upward at the rear part due to the buoyancy effect and the twin-vortex closure is successfully predicted as shown in Fig. 4(b), which is also observed in experiments.

Second, the numerical approach is validated by the ventilated supercavity in the unsteady flow settings. The same conditions with the steady validation adding a periodical flow set as $fg=10\u2009Hz$ and $\gamma g=4\u2009deg$ are employed in both the numerical simulation and the experiment. The numerical results (computational fluid dynamics, for short) and experimental data (EXP, for short) are shown in Fig. 5(a) by comparing the $Dmax$ and its locations ($L1/2$,$yDmax$). The location of maximum diameter is defined in the standard *yz*-coordinate system with the origin at the cavitator center. All parameters are normalized by the gust wavelength $\lambda g=U\u221e/fg$. The $L1/2$ and $yDmax$ in the experiments are measured through an automatic image processing technique described by Lee et al. [13] with an uncertainty around 8%. It is depicted that $Dmax$ and its location are fluctuated periodically both in simulation and in experiments. The simulated $Dmax$ and $yDmax$ vary synchronously with the experiments. In the simulation, the $Dmax$ fluctuates only about 1% around its mean, while the $L1/2$ fluctuation is about 16.5%, matching closely with the experiment values (i.e., 1.1% for $Dmax$, and 17.0% for $L1/2$, respectively). However, the fluctuation of $yDmax$ from the experiment is about 30% higher than that from the simulation. We attribute such discrepancy to the differences in the unsteady flow setup between the simulation and experiment. Specifically, the varying blockage associated with the hydrofoil pitching in the experiments leads to the stronger shape fluctuation of the cavity.

The supercavity pressure fluctuation in a gust cycle is also compared across the experiment and the simulation (Fig. 5(b)). The test section and cavity pressures captured in the simulation and the experiment are normalized by the dynamic pressure of the water flow. The experiment and simulation show similar periodical variation of cavity pressure and cavitation number. Additionally, both the experiment and the simulation show double peaks of the test section and cavity pressure variations within a gust cycle, corresponding to two instances of the hydrofoils with the maximum flapping angles. It is worth noting that the amplitude of the fluctuation in the normalized test section pressure is around 0.2% which is larger than the uncertainty of the pressure measurement around 0.1% corresponding to 100 Pa. In the simulation, $p\u221e$ fluctuates about 0.5% around its average, and $pc$ has a fluctuation around 0.6%, in comparison with their counterparts of the experiments results both around 0.2%. Nevertheless, the fluctuation of $\sigma c$ for the simulation is around 0.3% around its average value and showing a slight phase difference with the experimental result fluctuating at the amplitude of 2.4%. The varying minor loss of the tunnel corresponding to different phases of the flapping foils attenuates the amplitudes of pressure fluctuations during the experiments. This minor loss also leads to the varying tunnel velocity and a local minimum of $p\u221e$ at the minimum angle of attack of the foils and a local maximum corresponding to the maximum angle of attack of the foils (note that the tunnel is vented to the atmosphere during the experiments). Therefore, an out-of-phase response of the $p\u221e\u2009$ and $pc\u2009$ to the phase change of flapping foils occurs and eventually causes the enhanced $\sigma c$ fluctuations during the experiments.

## 3 Results

Based on the validations in Sec. 2.3, different unsteady flow conditions listed in Table 2 are simulated with periodic gust flows corresponding to different gust generator setups ($fg$ and $\gamma g$) at fixed $\u2009dc=20\u2009mm$, $Fr=18.7,\u2009\u2009CQ=0.15$, and $dc=20\u2009mm$. The transient supercavity behavior and the cavity shape deformation are provided during one cycle for the case of AOA8f10 state ($fg=10\u2009Hz$ and $\gamma g=8\u2009deg$), and the effects of $\gamma g$ and $fg$ are discussed subsequently. For all the simulations, the cavity boundaries are defined as the isosurface of the air volume fraction of 10% ($\alpha a=0.1$). It is noted that quantification of the vertical velocity component and the wavelength nondimensionalized by the cavity dimension are listed in Table 2. These parameters are correlated with different sea states [13]. As shown in Table 2, our work could cover sea states 2, 3, 4, which is corresponding to >50% of the possible sea states according to Ref. [13]. The sea states are used to characterize the influence of surface wave-induced fluctuation on the vehicles operating below sea surfaces. Our study provides a fundamental understanding of such influence, which is generally applicable to vehicles operating near the free surfaces including surface shape with lubrication (i.e., not limited to high-speed torpedoes).

State | Angle of attack (deg) | $fg$ (Hz) | $vgmax$ (m/s) | λ (m) | Sea state [13] |
---|---|---|---|---|---|

AOA8f10 | 8 | 10 | 0.227 | 0.85 | 4 |

AOA6f10 | 6 | 10 | 0.327 | 0.85 | 3 |

AOA4f10 | 4 | 10 | 0.448 | 0.85 | 3 |

AOA8f5 | 8 | 5 | 0.465 | 1.70 | 3 |

AOA8f1 | 8 | 1 | 0.5 | 8.50 | 2 |

State | Angle of attack (deg) | $fg$ (Hz) | $vgmax$ (m/s) | λ (m) | Sea state [13] |
---|---|---|---|---|---|

AOA8f10 | 8 | 10 | 0.227 | 0.85 | 4 |

AOA6f10 | 6 | 10 | 0.327 | 0.85 | 3 |

AOA4f10 | 4 | 10 | 0.448 | 0.85 | 3 |

AOA8f5 | 8 | 5 | 0.465 | 1.70 | 3 |

AOA8f1 | 8 | 1 | 0.5 | 8.50 | 2 |

### 3.1 Transient Supercavity Behaviors.

Figure 6 shows the supercavity variation at eight instants during a gust cycle for AOA8f10 state ($fg=10\u2009Hz$ and $\gamma g=8\u2009deg$). Inset figures of $\alpha a$ overlaid with the velocity vector field near the closure region are also provided at selected instants. At $t=0\u2009T$, where $T$ is the period of the foil flapping, the supercavity encompasses the cavitator, grows along the strut, and closes with two hollow vortex tubes. The cavity body shows an upward shift at the rear region due to the buoyancy effect. In accordance with periodic flow unsteadiness, the supercavity body fluctuates periodically in the vertical direction at the symmetry plane, especially at its closure region. In addition, it exhibits a periodical growth/shrink in the streamwise direction with the cavity length illustrated by the dash lines. Specifically, during $t=0\u223c2/8\u2009T$, the cavity body extends downstream, and meanwhile the two hollow vortex tubes moves upward and enlarges. Both the cavity length and the diameter of vortex tubes reach the maximum values at $t=2/8\u2009T$. Subsequently, the cavity body begins to shrink along the ventilation pipe until the minimum length at $t=6/8\u2009T$. As a result, the twin vortex tubes reduce in diameter and then disappear at $t=6/8\u2009T$, transitioning to a reentrant jet closure. Finally, both the cavity and twin vortex tubes grow again at $t=7/8\u2009T$ followed by a new cycle. Moreover, a low-$\alpha a$ portion is clearly observed near the cavity closure at $t=0\u2009T$ and $t=1/8\u2009T$, indicating an entrant-jet type of flow. Such flow declines during $=2/8\u223c6/8\u2009T$ and reoccurs at $t=7/8\u2009T$. As shown in the close-up view in Fig. 6 at $t=2/8\u2009T$, near the cavity closure region, the majority of the flow field within the vertical extent of the cavity is occupied by a reverse flow except in the internal boundary layer, which refers to a flow region of high shear extended from the liquid–gas interface within the body of the cavity. Wu et al. [11] observed such internal boundary layers in the upstream portion of a supercavity generated by a backward-facing cavitator (i.e., the cavitator with the mounting strut in front of the cavitator) and postulated the similar internal flow field near the closure region. Additionally, in comparison with Wu et al. [11], the simulation uses different cavitator configuration with a ventilation pipe behind the cavitator. As a result, as depicted in the close-up views at $t=2/8\u2009T$ and $t=4/8\u2009T$, a velocity boundary layer develops on the surface of the ventilation pipe with strong reverse flow.

The variation of internal flow field along the supercavity is further examined with distribution of air volume fraction and vorticity at cross section plane at different streamwise locations as shown in Fig. 7. The distributions of $\alpha a$ and the streamwise vorticity component ($\Omega z=\u2202uy/\u2202x\u2212\u2202ux/\u2202y$) are illustrated in four cross sections (termed as S1, S2, S3, and S4). Note that the S1 is the location of the cavity maximum diameter and others are near the closure region. The supercavity is symmetric at the S1 section and so is the internal flow depicted by the vorticity distribution. At S2 and S3 sections, the buoyancy effect causes the cavity to curve upward and its cross section develops to a “crescent” shape with high volume fraction as reported in Ref. [18], which generates two pairs of counter-rotating vortices at S2 and S3 sections. At the section S4, the air volume fraction decreases dramatically, and the outer pair of vortices becomes dominant with the inner pair of vortices located around the ventilation pipe. Such vortex pairs are associated with the velocity field at the cross section as depicted in Fig. 8, which shows the air contours corresponding to the air volume fraction about 10% along with the vortex outlines. The outer vortex pair is generated from the strong velocity gradient near the water–gas interface of the cavity, while the inner pair of vortices is a result of the velocity boundary layer developed along the mounting strut, which is not present in the internal flow experiment from Ref. [11] due to different cavitator configurations.

### 3.2 Cavity Shape Deformation.

First, the variation of the cavity overall geometry is depicted as Fig. 9 by showing the time-dependent $Dmax$ of the supercavity and its locations ($L1/2,\u2009yDmax$). The cavity half length $L1/2$ varies at twice of $fg$ with the peaks located at the extrema of the maximum hydrofoil vertical speeds, while the frequency of $Dmax$ and $yDmax$ stays the same as that of the gust flow. Additionally, $L1/2$ oscillates at an amplitude around 33.0% comparing to its average while $Dmax$ is around 3.5%. Comparatively, $yDmax$ has a fluctuating amplitude around 64.3% relative to its maximum value. The variation of cavity morphology is characterized by the change of its instantaneous centerline at the symmetric *yz*-plane extracted from its isosurface at eight consecutive time instances ($t=0\u223c7/8\u2009T$) in a gust cycle. It is worth noting that the cavity is located in the wavy state described in the previous experimental investigation [15]. The streamwise position *z* is normalized by the gust wavelength (i.e.,$\u2009\lambda g=U\u221e/fg=0.85\u2009m$), and vertical location is normalized by the cavitator size ($\u2009dc$) in Fig. 10. The cavity centerline exhibits wavelike deformation in the unsteady flow similar to the experimental results. It is worth noting that due to the buoyancy, the cavity has relatively stronger fluctuation in its centerline after $L1/2$_{.} The minimum vertical position of the centerline is advected at the speed of 96.5% of the freestream flow speed. Additionally, the amplitude of fluctuation of the cavity centerline in *y*-direction comparing to its neutral location is around 6 mm or 0.7% relative to the wavelength which is in the same order of the wave amplitude (0.8% of relative to the wavelength). Similar cavity geometrical results (i.e., cavity deformation advection speed and degree of surface deformation) are observed from the experiments under same flow conditions [15]. Therefore, the numerical settings presented in this paper are able to capture the geometrical behaviors of wavy state cavity in the unsteady flow.

### 3.3 Pressure Distributions and Ventilated Cavitation Number.

The instantaneous cavitation number is calculated from the difference between test section pressure and cavity pressure normalized by the dynamic pressure of the incoming flow as shown in Fig. 11. Compared to their averaged values, dimensionless test section pressure ($P\u221e$) has a fluctuation of 1.9% around its mean, and $Pc$ fluctuates slightly stronger at 2.2%. The cavity pressure having a phase delay of 0.014 s with respect to the phase of hydrofoils is associated with the propagation of the perturbation induced by the flapping hydrofoils. Figure 12 shows the pressure coefficient, $Cp=2(p\u2009\u2212p\u221e)/(\rho wU\u221e2)$, calculated at eight instants during a gust cycle at AOA8f10 state, to give an insight into the pressure distribution inside the cavity. Similar to the steady case [17], in a gust cycle, the pressure distribution remains largely uniform in the majority portion of the cavity at the front with a variation of pressure gradient near the cavity closure. Near the closure region, the pressure shows a steeper gradient along the streamwise direction at $t=1/8\u2009T$ and $2/8\u2009T$ compared to other phases in a cycle, as the low-pressure air inside the cavity rapidly increases its pressure to match the high pressure of the water after being discharged from the hollow vortex tube [24]. Such change of pressure gradient at the closure leads to the closure modes variation as reported in Refs. [14] and [15]. Specifically, as the closure mode transitioning from reentrant jet closure (corresponding to $t=6/8\u2009T$) to twin vortex closure ($t=2/8\u2009T$), the pressure at the closure of the cavity first decreases sharply and then increases as the air moves downstream into the vortex tube. The connection between the pressure variation and the change of closure modes during the gust cycle can be explained using the framework proposed in Karn et al. [10] by introducing a differential pressure term, $\Delta P\u0303=(pout\u2212pin)/(1/2\rho wU\u221e2)=P\u0303out\u2212P\u0303in$, where $P\u0303in$ and $P\u0303out$ represent the dimensionless static pressure inside and just outside the cavity at the closure, respectively. The $P\u0303in$ is calculated to be the average along a vertical line corresponding to a mean air volume fraction of 0.9 in the vicinity of cavity closure location (marked by the black dashed line in Fig. 6), whereas $P\u0303out$ is the averaged pressure in the supercavity wake. As shown in Table 3, the $\Delta P\u0303$ is 0.36 for a reentrant jet type closure at $t=6/8\u2009T,$ which is considerably larger than the corresponding value of 0.12 for a twin vortex closure at $t=2/8\u2009T$, providing further support to the pressure criterion proposed by Karn et al. [10].

Time | $P\u0303in$ | $P\u0303out$ | $\Delta P\u0303$ |
---|---|---|---|

0/8T | 0.20 | 1.87 | 0.20 |

1/8T | 0.15 | 1.86 | 0.13 |

2/8T | 0.09 | 1.85 | 0.12 |

3/8T | 0.17 | 1.90 | 0.17 |

4/8T | 0.19 | 1.91 | 0.19 |

5/8T | 0.28 | 1.96 | 0.28 |

6/8T | 0.36 | 2.01 | 0.36 |

7/8T | 0.21 | 1.88 | 0.21 |

Time | $P\u0303in$ | $P\u0303out$ | $\Delta P\u0303$ |
---|---|---|---|

0/8T | 0.20 | 1.87 | 0.20 |

1/8T | 0.15 | 1.86 | 0.13 |

2/8T | 0.09 | 1.85 | 0.12 |

3/8T | 0.17 | 1.90 | 0.17 |

4/8T | 0.19 | 1.91 | 0.19 |

5/8T | 0.28 | 1.96 | 0.28 |

6/8T | 0.36 | 2.01 | 0.36 |

7/8T | 0.21 | 1.88 | 0.21 |

### 3.4 Effect of the Angle of Attack.

The effect of the angle of attack on cavity behavior is investigated through the simulation under three angles of attack, i.e., $\gamma g=4,\u20096,\u20098\u2009deg$, with a fixed gust frequency at $fg=10\u2009Hz$ at the same inlet flow speed and ventilation rate. Figure 13 shows the deformation of the cavity centerline in a gust cycle upon different angles of attack. Compared the results shown in Fig. 10 (the result is also added in Fig. 13 for better comparison), under lower angles of attack, the cavity still exhibits wavelike undulation. However, with the reduction of $\gamma g$, there is a corresponding decrease of the undulation amplitude along the entire cavity. At $\gamma g=8\u2009deg$, the cavity centerline has the maximum undulation at $t=0\u2009T$ at the end of the cavity (i.e., corresponding to $z/\lambda g=0.4$ in the figure). Remarkably, as the $\gamma g$ changes from $8$ to $4\u2009deg$, the undulation amplitude of the cavity centerline at $t=6/8\u2009T$ gradually catches up and then exceeds that at $t=0\u2009T$ at the end of the cavity. Figure 14 presents the variation of test section pressure and cavity pressure in a gust cycle under the three angles of attack. As it shows, the pressure signals still exhibit periodic variation and a double-peak pattern in a gust cycle, but with a reduced amplitude corresponding to the decrease in the angle of attack. It is worth noting that there is a slight shift of the undulating phase of the pressure signals across different cases of $\gamma g$ (i.e., instantaneous pressure of high angle of attack displays larger delaying to the phase change of the hydrofoils comparing to smaller angles of attack cases). We attribute such shift to the increasing blockage effect due to the stronger shape deformation for high angle of attack case. In addition to the cavity shape and the pressure signals, the internal flow of the supercavity exhibits similar features at different cases of $\gamma g$ with increasing reverse flow amplitude and increasing thickness of the boundary layers for high $\gamma g$ (not shown for brevity).

### 3.5 Effect of the Gust Frequency.

The effect of gust frequency is studied through the simulation under three gust frequencies, i.e., $fg=1\u2009Hz,\u20095\u2009Hz,\u200910\u2009Hz$, with a fixed angle of attack at $\gamma g=8\u2009deg$ at the same inlet flow speed and ventilation rate. With increasing $fg$, the gust wavelength $\lambda g=U\u221e/fg$ and wave amplitude characterized by maximum vertical velocity of the disturbance ($vgmax$) decrease. When the gust wavelength becomes comparable to the cavity length (corresponding to $fg=10\u2009Hz$), the overall cavity centerline exhibits clear wavelike undulation (Fig. 15). Such undulation behavior diminishes significantly with reducing $fg$. For all three $fg$, the pressure signals fluctuate periodically with a double-peak in a gust cycle (Fig. 16). However, both test section pressure and cavity pressure have stronger fluctuations with increasing $fg$ associated with the augmented blockage effect and shape deformation at higher $fg$. Similar to the case of increasing $\gamma g$, the pressure signals for $fg=\u200910\u2009Hz$ shows a slightly larger delay to hydrofoil phase in comparison with the case of 1 Hz. This phase difference is likely to be caused by the larger cavity deformation for $fg=\u200910\u2009Hz$ causing increased blockage that lowers the advection speed of the incoming flow perturbation in a closed tunnel. In contrast, the cavitation number at lower $fg$ displays much stronger fluctuation. Such trend can be explained from the comparison of internal flow results across three cases. Under the cases of lowering $fg$, the cavity has remarkably stronger reverse flow in amplitude and a thicker boundary layer. Correspondingly, although the cavity shape variation diminishes with reducing $fg$, the cavity exhibits a more prominent closure variation associated with larger cavitation number fluctuation. Such trend of cavity closure variation under different $fg$ has also been reported in experiment [15] and can be explained by the theoretical analysis in Refs. [10] and [14] which connect dimensionless differential pressure $\Delta P\u0303$ introduced in Sec. 3.33B to the cavitation number.

It is worth noting that the phase of pressure fluctuation for $fg=5\u2009Hz$ is ahead of the other two cases. We suggest that such difference is related to the phase difference of cavity shape variation under different gust frequencies. Specifically, for $fg=5\u2009Hz$, the instant of maximum cavity shape deformation is at $t=6/8\u2009T$ while that for $fg=10\u2009Hz$ takes place at $t=7/8\u2009T$ (Fig. 15). Therefore, for $fg=5\u2009Hz$, the influence of blockage kicks in ahead that for $fg=10\u2009Hz$, leading to the phase difference in their pressure signals. However, for $fg=1\u2009Hz$, the overall shape variation of cavity is negligible comparing to the other cases and thus does not alter the phases of the pressure. In comparison with the experimental results under the same flow conditions and cavitator model reported in Ref. [15], the simulation does not capture the change of supercavity mode as it transitions from wavelike undulation to cavity pulsating with reduction of $fg$. Specifically, the experiment shows drastic change of cavity lengths due to the shed-off of gas pocket at the rear of the cavity. Such instability has not been captured in our simulation due to the limitation of the numerical approach. This limitation is discussed in detail in Sec. 4. Nevertheless, the simulation can still provide an overall degree and phase of wavelike undulation occurring on cavity surface under different unsteady conditions.

## 4 Summary and Discussions

The ventilated supercavitation from a forward facing cavitator under periodic gust flow is numerically studied by using the mixture model in the scope of the Eulerian–Eulerian framework. The flow unsteadiness is simulated by adding a harmonic varying vertical velocity term matching that produced from a pair of flapping hydrofoils in the experiments [13,15]. The numerical method is first validated through the experiments under specific steady and unsteady flow conditions. It has been shown that the simulation can capture the degree of cavity shape fluctuation and internal pressure variation in a gust cycle. Specifically, the cavity centerline shows periodic wavelike undulation with a maximum amplitude matching that of the incoming flow perturbation. The cavity internal pressure also fluctuates periodically, causing the corresponding change of difference between internal and external pressure across the closure that leads to the closure mode change in a gust cycle. In addition, the simulation captures the variation of cavity internal flow, particularly the development internal flow boundary layer along the cavitator mounting strut, upon the incoming flow perturbation, correlating with cavity deformation and closure mode variation. With increasing angle of attack, the cavity exhibits augmented wavelike undulation and pressure fluctuation. As the wavelength of the flow perturbation approaches the cavity length with increasing gust frequency, the cavity experiences stronger wavelike undulation and internal pressure fluctuation but reduced cavitation number variation.

The numerical simulation provides detailed information of cavity internal flow and pressure distribution which are very difficult to be measured even with the most recent experimental investigation employing sophisticated imaging techniques [11]. Specifically, the simulation result captures the development of internal flow throughout the entire cavity, including the region near the closure region which is not accessible in Ref. [11]. Further, in comparison with backward-facing cavitator used in Ref. [11], the simulation is conducted on a forward-facing cavitator and shows the change of internal flow structures due to the influence of the mounting strut, such as the enhanced reverse flow due to the formation of the boundary layer along the mounting strut. The internal pressure distribution provided from the simulation allows us to directly assess the variation of pressure difference across the closure region of the cavity. The result not only provides support to the framework proposed by Karn et al. [10] for explaining the supercavity closure change under different steady flow conditions, but also suggests the applicability of such framework for predicting cavity closure variation in unsteady flows.

Based on the simulation under different angles of attack and gust frequencies, we have shown that the phases of cavity deformation, particularly the location and the time instant corresponding to the maximum vertical displacement of the cavity, vary according to the wavelength and wave amplitude of unsteady flow. Such variation needs to be considered in the ventilation design and the control of the supercavitating devices. Additionally, the simulation shows that the varying degree of cavity deformation can cause phase shift of pressure signals under different unsteady conditions in a closed water tunnel. This observation further underscores the importance of understanding the effect from flow facilities on supercavity behavior as discussed in Ref. [27].

It is noteworthy that the mixture model is used for the numerical simulation which is in the Eulerian framework, but it is strictly a “one-fluid” model. For the present numerical approach, the multiphase flows are inhomogeneous, that is, there is a relative velocity between the water and the air, which is shown in Fig. 17. The relative velocity (*u*_{w}−*u*_{a}) distributions are depicted on S1, S2, and S3 cross sections of the cavity at $t=0T$ for the AOA8f10 state ($\gamma g=8\u2009deg$ and $fg=10\u2009Hz$) overlaid with air fraction contour corresponding to air volume fraction of 10% (marked in black). Gradients of the relative velocity are observed within the cavity, but the inner boundary layer is not resolved by the present simulations. Besides, the present simulation is unable to capture the cavity length change due to shed-off of gas pocket near its closure region (i.e., pulsation of the cavity) which occurs under the cases of AOA8f1 and AOA8f5 in the experiments [15]. Therefore, the Eulerian model (“two-fluid model”) is considered in our future work.

As discussed in Ref. [15], the transition from the wavelike undulation to the cavity pulsation is a result of the bubble breakup near cavity closure induced by the instability of the gas–liquid interface. Therefore, to fully capture cavity behaviors like pulsation, future numerical simulation should consider detailed bubble–liquid interaction such as bubble induced turbulence, and lift force and drift velocity of bubbles which are not explicitly included in the present simulation. As an example, in the modeling of a bubble column, Selma et al. [28] considered additional terms like lift force on individual bubbles, bubble drifting due to turbulence, and virtual mass force in the exchange of momentum on gas–liquid interface. A standard $k\u2212\epsilon $ model with bubble induced turbulence was used for turbulence modeling. The simulation results have shown to be able to resolve the vortex structures developed on the gas–liquid interface due to flow instability. Large eddy simulation with locally adaptive eddy viscosity terms provides another approach to capture the cavity pulsation in unsteady flows. It was recently employed to capture the pulsation of a hydrofoil cavitation originated from flow instability on the vapor–water interface [29]. Such high-fidelity simulation should be applied to the study of ventilated supercavitation in unsteady flows, which can provide more physical insights into its behaviors.

## Acknowledgment

Renfang Huang is supported by the State Scholarship Fund from China Scholarship Council and the National Key R&D Program of China (2017YFC1404200). We also acknowledge the support from the Office of Naval Research (Program Manager, Dr. Thomas Fu) under Grant No. N000141612755 and the start-up funding received by Professor Jiarong Hong from University of Minnesota.

## Funding Data

China Scholarship Council (Funder ID: 10.13039/501100004543).

National Key R&D Program of China (Grant No. 2017YFC1404200; Funder ID:10.13039/501100002855).

Office of Naval Research (Grant No. N000141612755; Funder ID: 10.13039/100000006).

University of Minnesota (Funder ID: 10.13039/100007249).

## Nomenclature

- $Am$ =
interfacial area per unit volume

- $CD$ =
nondimensional drag coefficient =0.44

*C*_{p}=pressure coefficient $=2(p\u2009\u2212p\u221e)/(\rho wU\u221e2)$

*C*_{Q}=air entrainment coefficient

*d*_{c}=cavitator diameter

- $Dmax$ =
maximum diameter of the supercavity

- $fg$ =
frequency of periodic flows (wavelength)

- Fr =
Froude number

*g*=gravitational acceleration

*L*_{1/2}=half-length of the supercavity

- $Mk$ =
total interfacial drag forces

*p*=pressure

*p*_{c}=cavity pressure

*p*_{∞}=ambient pressure upstream of the cavitator

- $P\u0303c$ =
dimensionless cavity pressure

- $P\u0303in,\u2009P\u0303out$ =
dimensionless static pressure inside and just outside the cavity at the closure

- $P\u0303\u221e$ =
dimensionless test section pressure

- $Q\u02d9$ =
volumetric ventilation rate

- $Sk,buoy$ =
momentum source due to buoyancy forces

- $uk$ =
mean velocity vector

*U*_{∞}=upstream incoming flow velocity

- $\alpha k$ =
volume fraction of phase

*k* - $\gamma g$ =
angle of attack of two hydrofoils

- $\Delta P\u0303$ =
dimensionless pressure difference

- $\lambda g$ =
gust wavelength, $\lambda g=U\u221e/fg$

- $\mu k$ =
molecular dynamic viscosity of phase $k$

- $\mu t,\u2009k$ =
turbulent eddy viscosity of phase

*k* - $\rho a$ =
air density =1.185 kg/m

^{3} - $\rho k$ =
density of phase

*k* - $\rho m$ =
density of the mixture

*ρ*_{w}=water density = 997 kg/m

^{3}*σ*=ventilated cavitation number

- $\u2009\u2009vg(t)$ =
vertical velocity term, $vg(t)=vgmax\u2009sin\u2009(2\pi fgt)$

- $vgmax$ =
maximum vertical velocity (amplitude)

- $\Omega z$ =
vorticity component

## References

*k*–

*ε*Model of Turbulence