**1. Introduction**

Cavitation is a complex phenomenon that occurs when the local pressure in a liquid drops below the local saturated vapor pressure. This pressure drop always occurs around a rapidly moving object and is generated at the internal or liquid-solid interface of the liquid [1,2]. Cavitation can cause performance degradation, corrosion, vibration, and noise in hydraulic machinery, and is a major factor in fluid instability [3–5]. Particularly for the underwater region of a high-speed vehicle located near a free surface, free surface-cavity interactions render the flow field more complicated [6].

Unsteady cavitation flows are an important problem for ocean engineering [7–9]. Knapp [10] observed and analyzed cloud-cavitation flows, demonstrating that the re-entrant flow starting from the closed position of the cavitation causes the cloud cavitation to break off. Kubota et al. [11] measured the flow structure of cloud cavitation using a laser Doppler anemometer (LDA) and compared it with the results obtained by a high-speed camera. The authors found a maximum of vorticity in the center of the cloud cavitation. Kawanami et al. [12] studied the nature of cavitation flow instability and found that cloud-like cavitations are closely related to retroreflective flows in a closed bubble. Reissman et al. [13] studied cloud-like cavitation around an oscillating hydrofoil using high-speed cameras and pressure sensors, and investigated the mechanism of transient pressure pulses near the detached cavitation. Arndt et al. [14] studied the flow of an NACA0015 hydrofoil from flaky cavitation to cloud cavitation by combining experimental and numerical methods. The authors observed two competing cavitation shedding induction mechanisms: re-entrant flow and a shock wave with the cavitating flow; moreover, they presented the dominant conditions for these two mechanisms. Gopalan et al. [15] proposed that the vortex generation and collapse are important characteristics of a cavitation flow. Leroux et al. [16,17] observed the cavitation flow of a NACA66 hydrofoil and proposed that the re-entrant flow is the primary cause of cavitation shedding. Passandideh-Fard and Roohi [18] used a modified VOF method to capture the gas-liquid interface of a two-dimensional cavitation flow field. The results show that the modified VOF method has higher calculation efficiency and calculation accuracy. Ganesh et al. [19] used wedges to study the transformation from sheet cavitation to cloud cavitation. Pendar and Roohi [20] used large eddy simulation to study cavitation around a sphere and compared with non-cavitation flow field. The result shows that cavitation suppresses instability near wake region and delays 3D decomposition of the vortex. Capurso et al. [21] simulated the hydrofoil with a passive control system. Their study indicated that passive control system suppress the cavitation effectively, which provides an effective solution to avoid the influence of cavity erosion on hydrofoil. Cheng et al. [22,23] used large-eddy simulations (LES) to capture an unstable tip leakage cavitation flow, and the spatial evolution of the tip leakage cavitation flow was divided into three stages. The results revealed how cavitation affects vorticity and turbulence. Researchers have performed numerous studies on cavitating flows at an infinite water depth; however, studies on cavity dynamics near free surface are limited.

Dawson [24] observed supercavitation near a free surface; however, the cavitating flow in the experiments was generated by ventilation under low-velocity restriction in a tunnel. Theoretical analysis and numerical simulations are the primary research methods used at present [25–28]. Faltinsen [29] proposed a theoretical model for free surface effects of hydrofoil cavitation. Kinzel et al. [30] simulated the supercavitation of the hydrofoil under a free surface and proved that using time accurate multiphase Navier-Stokes CFD method to simulate cavitation near the free surface is feasible. Wang et al. [31] performed experiments and simulations on the cavitation of a revolving body near a free surface. It has been proposed that the free surface will reduce the re-evaporation of the cavitation, and that the cavitation will become more stable as it approaches the free surface. Based on Wang's result, a multiphase flow solver for simulating near-free-surface cavitation flows has been proposed and validated [32,33]. The numerical simulation has been proven to be capable of analyzing the effect of unsteady cavitation flows near a free surface. However, the cavitation flow field evolution and cavitation dynamics under near-free-surface condition remain unclear.

Due to the limitations of experimental measurement techniques, increasingly more scholars have begun to use numerical methods to study cavitation flow phenomena. Due to the strong instability of the flow near the cavitation, the simulation method plays an important role in the numerical calculation. Most previous simulations used the Reynolds average Navier-Stokes (RANS) method for cavitation calculations due to its low computational complexity and high computational efficiency [34,35]. However, this method overestimates the turbulent eddy viscosity and the peak pressure pulsation of the cavitation [36]. With recent advances in computational power, LESs have arisen as a superior choice for simulating cavitation flows [37]. Pendar and Roohi [38] studied the cavities length and diameter of hemispherical heads and cone cavitators under different cavitation numbers by using different turbulence models and cavitation models. As a result, the combination of different turbulence models and cavitation models on the surface has different predictions of cavitation

shapes. Sedlar et al., [39] used an NACA2412 hydrofoil to simulate cavitation flows under different turbulence models. The results showed that the pressure pulsation and shedding cycle frequency obtained by the LES are similar to the experimental values. The LES model can predict the unstable nature of the flow, accurately capture small-scale vortices, forecast the flow field, and provides insight into the flow details [40–42]. Therefore, the cavitation-flow LES model has been used to predict multiscale turbulent vortex characteristics.

For complex cavitation-vortex interactions in cavitating flow, it is important to study the coherent structure of turbulence. In recent years, modal decomposition techniques have been gradually applied to computational fluid dynamics. These techniques are considered as important analytical tools for assessing nonlinear fluctuation data. Intrinsic mode decomposition (POD) and dynamic mode decomposition (DMD) are the two primary methods of modal decomposition. However, the POD method may contain multiple frequency components in single modality; thus, the frequency information of each coherent mode cannot be provided [43,44]. The DMD method sorts the components according to the frequency, with single-frequency modalities that show strong advantages in the analysis of flow characteristics. The standard DMD method is based on the Koopman analysis of nonlinear dynamical systems, which was originally proposed by Schmid et al. [45]. Vinha et al. [46] used the DMD method to analyze the jet pulsation behavior and cavitation flow caused by different jet shapes. Prothin et al. [47] and Liu et al. [5] combined POD with the DMD method to study the development of hydrofoil cavitation under a high Reynolds number. The cavitation flow field modes of two different positions of the hydrofoil were analyzed, and it was proposed that the DMD method can more accurately extract the frequency characteristics.

At present, most studies on unsteady cavitating flows over a hydrofoil have been performed for an infinite water depth. However, study on cavitation flows near the free surface is still limited, especially the mechanism of interaction between free surface and cavitation is unclear. This paper presents a numerical study of the unsteady cavitation dynamics over a NACA66 hydrofoil near the free surface. Section 2 introduces the numerical framework, including basic governing equations, cavitation models, turbulence models, numerical setup and validations. Section 3 presents the effects of the free surface on the unsteady cavitation. An improved DMD method is applied to investigate the unsteady cavitating flows fealtures. Finally, Section 4 summarizes the main conclusions.

#### **2. Numerical Methods**

#### *2.1. Governing Equations*

The volume fraction of phase *n* is:

$$
\alpha\_n = \frac{V\_n}{V} \tag{1}
$$

where *Vn* is the volume of phase *n* and *V* is the total volume.

The governing equations of the mixture model are:

$$\frac{\partial \rho}{\partial t} + \frac{\partial \{\rho u\_j\}}{\partial x\_j} = 0 \tag{2}$$

$$\frac{\partial(\rho u\_i)}{\partial t} + \frac{\partial(\rho u\_i u\_j)}{\partial x\_j} = -\frac{\partial p}{\partial x\_i} + \frac{\partial}{\partial x\_j} \left(\mu \frac{\partial u\_i}{\partial x\_j}\right) \tag{3}$$

where *p* is the pressure, and *ui* is the velocity along the *i*-th direction.

The mixture density ρ is:

$$
\rho = \alpha\_l \rho\_l + (1 - \alpha\_l)\rho\_\mathbf{v} \tag{4}
$$

where α*<sup>l</sup>* is the liquid volume fraction, ρ*<sup>l</sup>* is the liquid density and ρν is the vapor density.

#### *2.2. Large-Eddy Simulation Model*

The LES model can divide the vortices in a turbulent flow into large-scale vortices and small-scale vortices by using a filtering function. Direct numerical calculations are used to solve the large-scale vortices, and small-scale vortices are calculated via a connection with large-scale vortices. The LES governing equations are:

$$\frac{\partial \rho}{\partial t} + \frac{\partial (\rho \overline{u}\_j)}{\partial x\_j} = 0 \tag{5}$$

$$\frac{\partial(\rho \overline{u}\_i)}{\partial t} + \frac{\partial(\rho \overline{u}\_i \overline{u}\_j)}{\partial \mathbf{x}\_j} = -\frac{\partial \overline{p}}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \Big(\mu \frac{\partial \overline{u}\_i}{\partial \mathbf{x}\_j}\Big) - \frac{\partial(\rho \tau\_{ij})}{\partial \mathbf{x}\_j} \tag{6}$$

where τ*ij* is the nonlinear sub-grid scale (SGS) stress tensor, which is defined as:

$$
\pi\_{i\bar{j}} = \rho \overline{u\_i u\_{\bar{j}}} - \overline{u}\_i \overline{u}\_{\bar{j}} \tag{7}
$$

The product of the fluid strain rate, *Sij* and an assumed sub-grid viscosity υ*SGS* is used to describe the tensor, which is based on the Boussinesq hypothesis:

$$
\pi\_{i\dot{j}} - \frac{1}{3}\pi\_{kk}\delta\_{i\dot{j}} = -2\nu\_{\text{SGS}}\overline{\mathfrak{S}}\_{i\dot{j}}\tag{8}
$$

The SGS model used in this work is the WALE model [48], which is given as:

$$\nu\_{\rm SGS} = \rho L\_s^2 \frac{\left(S\_{ij}^d S\_{ij}^d\right)^{3/2}}{\left(\overline{S}\_{ij} \overline{S}\_{ij}\right)^{5/2} + \left(S\_{ij}^d S\_{ij}^d\right)^{5/4}} \tag{9}$$

$$\overline{S}\_{i\bar{j}} = \frac{1}{2} \left( \frac{\partial \overline{u}\_{\bar{i}}}{\partial \mathbf{x}\_{\bar{j}}} + \frac{\partial \overline{u}\_{\bar{j}}}{\partial \mathbf{x}\_{\bar{i}}} \right) \tag{10}$$

$$\delta S\_{ij}^d = \frac{1}{2} (\overline{\mathcal{g}}\_{ij}^2 + \overline{\mathcal{g}}\_{ji}^2) - \frac{1}{3} \delta\_{ij} \overline{\mathcal{g}}\_{kk'}^2 \overline{\mathcal{g}}\_{ij} = \frac{\partial \overline{u}\_i}{\partial \mathbf{x}\_j} \tag{11}$$

where *Ls* is the length or grid-filter width, which is given as:

$$L\_s = \min(kd\_\prime \mathcal{C}\_w V^{1/3}) \tag{12}$$

Here, *k* is von Karman's constant, *d* is the distance to the closest wall, *V* is the cell volume, and *Cw* is a model coefficient which is set to 0.5.

### *2.3. Cavitation Model*

Cavitation models generally consist of transport equations. These models are usually obtained by simplifying and improving the Rayleigh-Plesset cavity dynamics equation [49,50]. The processes of evaporation and condensation can be expressed by inserting a source term into the transport equation to describe the phase transition rate. The Schnerr-Sauer model is used as the cavitation model [49], and the transfer equation is:

$$\frac{\partial(\rho\_v \alpha\_v)}{\partial t} + \frac{\partial(\rho\_v \alpha\_v \mu\_j)}{\partial \mathbf{x}\_j} = \dot{m}^+ - \dot{m}^- \tag{13}$$

where *<sup>a</sup>*<sup>ν</sup> is the vapor volume fraction. The source terms . *<sup>m</sup>*<sup>+</sup> and . *m*<sup>−</sup> describe the evaporation and condensation which can be solved as:

$$\dot{m}^{+} = \frac{\rho\_{\upsilon}\rho\_{l}}{\rho} \alpha\_{\upsilon} (1 - \alpha\_{\upsilon}) \frac{3}{R\_{b}} \sqrt{\frac{2}{3} \frac{\max(p\_{\upsilon} - p\_{\prime}, 0)}{\rho\_{l}}} \tag{14}$$

$$\dot{m}^{-} = \frac{\rho\_{v}\rho\_{l}}{\rho} \alpha\_{v} (1 - \alpha\_{v}) \frac{3}{R\_{b}} \sqrt{\frac{2}{3} \frac{\max(p - p\_{v}, 0)}{\rho\_{l}}} \tag{15}$$

where . *<sup>m</sup>*<sup>+</sup> is the evaporation term, . *m*<sup>−</sup> is the condensation term, and *p*<sup>ν</sup> is the saturation vapor pressure at the local temperature. *Rb* is the cavity radius which can be solved as:

$$R\_b = \left(\frac{\alpha\_v}{1-\alpha\_v} \frac{3}{4\pi} \frac{1}{N\_b}\right)^{\frac{1}{3}} \tag{16}$$

where *Nb* is the density of cavity number. According to Schnerr and Sauer [51], *Nb* = 1013.

#### *2.4. Numerical Setup and Validation*

For all the numerical simulations in the present study the Star-CCM+ software was used. A NACA66 hydrofoil was used in this work, which has a chord length *C* = 0.15 m and a span of 0.192 m. The angle of attack of the hydrofoil is 8◦. The computational domain is shown in Figure 1. The velocity inlet of the computational domain is approximately 4*C* from the leading edge of the hydrofoil and the inflow velocity is *<sup>V</sup>*<sup>∞</sup> <sup>=</sup> 5.33 m/s, with a Reynolds number of *Re* <sup>=</sup> 0.8 <sup>×</sup> 106. The outlet pressure can control the reference pressure of the flow field to adjust the cavitation number. The cavitation number is defined as:

$$
\sigma = (p - p\_{\overline{v}}) / (0.5 p\_l V\_{\infty}^2) \tag{17}
$$

where *p* is the far field pressure, which is set by the outlet boundary conditions. In order to be consistent with the experimental conditions [16], σ is set to 1.25 in this work. The saturated vapor pressure of water is *pv* = 2367 Pa for a flow-field temperature of 20 ◦C. To capture the effect of cavitation on the tail of the free surface, the pressure outlet was set to 20*C* downstream of the leading edge.

**Figure 1.** Computational domain and boundary condition.

Figure 2a shows the mesh of the computational domain. The mesh is generated by the Star-CCM+ auto trimmed mesher and prism layer mesher. Transition area between prism layer mesh and trimmed mesh is unstructured type, and structured mesh is generated in other regions. The total number of mesh elements is approximately 3 <sup>×</sup> 10<sup>6</sup> and the y+ < 1 over the hydrofoil surface to meet quality criterion requirements for LES, as shown in Figure 2b.The time step for the calculation is set to <sup>Δ</sup>*<sup>t</sup>* = 2.5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> to ensure stability in the numerical calculation.

**Figure 2.** Mesh generation and y+ at *t* = 0.566 s around the hydrofoil. (**a**) Mesh generation at Mid-section and a visualization close to the profile, (**b**) y+ at *t* = 0.566 s.

Three pressure-monitoring points were set on the mid-section of the hydrofoil, as shown in Figure 3.

**Figure 3.** Locations of monitoring points.

The middle section along the spanwise length of the flow field is shown in Figure 4. The zero point of the primary coordinate system, XOY, is located at the leading edge of the hydrofoil. To represent the shape of the free surface, a local coordinate system, X1OY1, is established at the free surface of the velocity inlet.

**Figure 4.** Definition of the coordinate systems in the computational domain.

In this paper, six working conditions are set under different depth conditions including the absence of a free surface condition and *h* = 1.25*C*, 1*C*, 0.75*C*, 0.5*C*, and 0.25*C*. In addition to the difference in depth from the free surface, the reference pressure also changes with the depth of the hydrofoil to keep the cavitation number at a constant.

To validate the numerical model, results from the numerical simulation and a previous experiment [16] were compared. In this experiment, numerical videos were used to record the cavitating flow, numbers of piezo-resistive transducers were used to measured pressure. Figure 5 shows the cavitation shape for a single cavitation cycle. *T* represents the period of cavity shedding. The numerical simulation accurately captured the periodic characteristics of the hydrofoil cavitation, such as the growing, shedding and collapsing processes. The cavitation shape is also well predicted by the numerical simulation. It shows that the numerical cavitation shape matches well with the experimental cavitation shape.

**Figure 5.** Numerically predicted cavity evolution compared with experimental results [16].

A collapse of the cavitation can lead to dramatic fluctuations in pressure. Thus, the dynamic pressure is difficult to predict. Numerical predictions of the pressure fluctuations at points P1 and P2 are plotted in Figure 6. As indicated by the experimental measurements, the pressure fluctuations caused by the cavitation collapse and the periodicity of the cloud cavitation flow are well predicted. The numerical prediction of cavitation shape and pressure results are highly consistent with experimental results indicating that the numerical method can effectively predict the unsteady cavitation dynamics.

**Figure 6.** Pressure fluctuation on the suction surface of the hydrofoil at P1 and P2.

### **3. Results and Discussion**

#### *3.1. E*ff*ect of a Free Surface on the Dynamic Cavity Evolution and Hydrodynamic Load*

To evaluate the effect of a free surface on cavitating flows, Figure 7 shows the cavity shapes at eight typical times in one cycle for the cases of without free surface and near the free surface (*h* = 0.25*C*), where *T*<sup>0</sup> and *T*<sup>1</sup> represent the evolution cycles of without and with free surfaces, respectively.

It can be seen that the overall characteristics of cavitation evolution in a cycle in both cases are similar, including cavitation development (from *t*<sup>0</sup> to *t*<sup>0</sup> + 4*T*0/8), shedding (from *t*<sup>0</sup> + 5*T*0/8 to *t*<sup>0</sup> + 6*T*0/8) and collapsing (from *t*<sup>0</sup> + 6*T*0/8 to *t*<sup>0</sup> + 7*T*0/8). However, it is noted that the maximum length of the attached cavity for the near free surface condition is significantly shorter and the scale of cloud cavitation is smaller. Hence, the existence of free surfaces suppresses the cavitation intensity.

(**a**) **Figure 7.** *Cont.*

**Figure 7.** Comparison of the cavity evolution for a typical cycle for no free surface case and near the free surface case (*h* = 0.25*C*). (**a**) Cavity evolution in a typical cycle for the no free surface case, (**b**) Cavity evolution in a typical cycle for the near free surface case (*h* = 0.25*C*).

(**b**)

To explain the effect of free surface on the intensity of cavitation, we introduce the local cavitation numbers, which is defined as:

$$
\sigma\_{\text{local}} = (p\_{\text{local}} - p\_{\upsilon}) / (0.5 \rho\_l V\_{\infty}^2) \tag{18}
$$

where σlocal is the local cavitation number and *p*local is the local absolute pressure.

Figure 8a shows the local cavitation number and Figure 8b shows the local velocity field on the mid-section plane of the computational domain. The local cavitation number near the cavitation zone is close to zero. As the distance from cavitation increases, the local cavitation number increases, indicating that a cavitation is more difficult to produce. Changes in the cavitation number represent changes in absolute pressure. It can be seen from the corresponding local velocity field that the high speed area under no free surface condition is larger than near free surface condition. This indicates that the existence of the free surface reduces the velocity of the flow field around the hydrofoil and resulting in a weaker cavitation intensity. Two straight lines *L*<sup>1</sup> and *L*<sup>2</sup> are set along the y-direction at *x*/*C* = 0.2 and 0.3. The pressure value of *L*<sup>1</sup> and *L*<sup>2</sup> as the *x* coordinates, and the vertical distance from the hydrofoil head is the *y* coordinate, as shown in Figure 8b. When the value of *y* is small, the pressure is small due to the cavitation. The pressure increases as y increases. When *y* = 0.05 m, the pressure under the near-free surface condition is larger than no free surface condition. The results show that the pressure increases faster near free surface. Near the suction surface of the hydrofoil, the free surface makes the pressure change more drastic, which leads to a larger local cavitation number. Subsequently, it is difficult to generate a cavitation, and the cavity length is reduced.

**Figure 8.** Comparison of the local cavitation number (**a**) and local velocity field (**b**) on the mid-section plane without free surface and with the free surface, (**c**) vertical pressure distribution in hydrofoil cavitation region.

The effect of the free surface and the cavitation is also reflected in the pressure fluctuation on the hydrofoil surface, as shown in Figure 9a. As the cavitation develops, the pressure is equal to the saturated vapor pressure, and the pressure is stable. When the cavities begin to break off, the cloud cavitation produces a dramatic pressure fluctuation, which causes the pressure to increase.

This feature is reflected in the periodic fluctuations of the pressure curve. The cavitation shedding frequency for the near free surface condition is approximately twice that of the no free surface condition. Figure 9b,c show the power spectral density (PSD) of the pressure curve for the two cases. The PSD curve indicates that the dominant frequency of periodic cavitation shedding is approximately 17.35 Hz in the absence of a free surface and approximately 35.75 Hz near the free surface. This difference is due to the reduced cavity length that occurs near a free surface and the time required for the cavities to develop and collapse which leads to an increased cavitation cycle frequency.

**Figure 9.** Pressure fluctuations and the corresponding PSD at P3. (**a**) Pressure fluctuation on the suction surface of the hydrofoil at P3; (**b**) PSD of the pressure fluctuation at P3 for the no free surface case; (**c**) PSD of pressure fluctuation at P3 for the near free surface case.

Figure 10 shows fluctuations in the lift and drag coefficients, which are defined as:

$$C\_L = \frac{\text{Lift}}{0.5\rho \times V\_{\infty}^2 \times \text{C} \times \text{Span}} \tag{19}$$

$$\mathcal{C}\_{D} = \frac{\text{Drag}}{0.5\rho \times V\_{\infty}^{2} \times \text{C} \times \text{Span}} \tag{20}$$

**Figure 10.** Comparison of the lift and drag coefficient fluctuations for the no free surface case and the near free surface case. (**a**) Lift coefficient; (**b**) Drag coefficient.

Due to the periodic fluctuations in pressure, the fluctuations in the lift and drag coefficients over time also appear cyclical. However, due to the strong instability of the cavitation flow, the periodic fluctuations of the lift and drag coefficients are not obvious.

Comparison shows that both coefficients are smaller under the near free surface case than under the no free surface case, which indicates that the free surface weakens the lift and drag of the hydrofoil by reducing the cavity length. Both the lift and drag coefficients have the same influencing factors: free surface and cavitation. The lift and drag coefficients increase with the increase in submergence ratio *h*/*C* [29]. Another influencing factor is that the low-pressure region of the hydrofoil suction surface is also reduced due to the reduction in the length of the cavitation. As the length of the cavity decreases, the pressure difference between the suction and pressure plane of the hydrofoil decreases. Hence, the existence of free surfaces reduces the lift and drag coefficients compared with the case of no free surface.

#### *3.2. Analysis of the Flow Structure Based on the DMD Method*

When the cloud cavity breaks off, a strong vortex is generated. The analysis of vortex structures is an important component in studying cavitation flows. The DMD method can effectively extract the modal information of each order in the cavitation flow field, and can analyze the vortex structure of each mode. For details on the standard DMD method, we refer the reader to previous publications [45]. It is assumed that the total number of snapshots is *N* and the time-space flow field of an unsteady flow can be represented as a matrix *VN* **1** :

$$\mathbf{V\_1^N} = [\mathbf{v\_1}, \mathbf{v\_2}, \dots, \mathbf{v\_N}] \tag{21}$$

where the column vector *vi* represents the *i*-th snapshot of the flow field. Next, *A* is defined as an approximation of the nonlinear mapping, indicating the evolution of the matrix from one timepoint to the next, which is expressed by:

$$
\mathbf{v}\_{i\leftrightarrow 1} = A\mathbf{v}\_i \tag{22}
$$

Then, the dynamic flow system is represented by the mapping matrix *A*. Matrix *VN* **<sup>1</sup>** can be defined as:

$$\mathbf{V\_1^N} = \begin{bmatrix} v\_{1'} A v\_{1'} A^2 v\_{1'} \dots \dots \angle A^{N-1} v\_1 \end{bmatrix} \tag{23}$$

*J. Mar. Sci. Eng.* **2020**, *8*, 341

Assuming that *a*<sup>i</sup> are coefficients of *A*, the last snapshot *vN* can be formulated by the previous *N* − 1 vectors:

$$
\Delta \boldsymbol{\sigma}\_{\text{H}} = a\_1 \boldsymbol{\sigma}\_1 + a\_2 \boldsymbol{\sigma}\_2 + \dots + a\_{N-1} \boldsymbol{\sigma}\_{N-1} + r \tag{24}
$$

or:

$$\mathbf{A}\mathbf{V}\_1^{\mathbf{N}-1} = \mathbf{V}\_1^{\mathbf{N}-1}\mathbf{S} + \mathbf{r}\mathbf{e}\_{\mathbf{N}-1}^T \approx \mathbf{V}\_1^{\mathbf{N}-1}\mathbf{S} \tag{25}$$

where *r* is the residual vector, *e<sup>T</sup> <sup>N</sup>*−<sup>1</sup> is the (*<sup>N</sup>* <sup>−</sup> 1)-th unit vector and matrix *<sup>S</sup>* is formed as:

$$\mathbf{S} = \begin{bmatrix} 0 & & & a\_1 \\ 1 & 0 & & & a\_2 \\ & & \ddots & \ddots & \vdots \\ & & & 1 & 0 & a\_{N-2} \\ & & & & 1 & a\_{N-1} \end{bmatrix} \tag{26}$$

when the residual *r* is small, the eigenvalue of *S* approximates the eigenvalue of *A*. The *S* matrix can be regarded as a low-dimensional form of matrix *A*; therefore, the eigenvalue of *S* can represent the main eigenvalues in matrix *A*. Eigenvalue decomposition is performed on matrix *S*:

$$\mathbf{S} = \mathbf{P}\mathbf{N}\mathbf{P}^{-1}, \mathbf{N} = \text{diag}(\mu\_1, \dots, \mu\_{N-1}) \tag{27}$$

where *P* is the matrix of the eigenvectors of matrix *S*.

The dynamic mode ϕ can be formed as follows:

$$
\mathfrak{sp} = \mathsf{V}\_1^{\mathsf{N}-1} A \tag{28}
$$

The first *m* flow field snapshots can represent the flow field snapshot at any time:

$$w\_{l} = \sum\_{j=1}^{m} \mu\_{l}^{j-1} q\_{j} \tag{29}$$

where *m* is the number of DMD modes and ϕ*j* is the column vector of the matrix ϕ.

The standard DMD method is not suitable for the convergence of strong unsteady flows such as cavitation [43]. Thus, in this paper, the raw velocity data are first processed, and then, DMD decomposition is performed.

The Hilbert-Huang transform is performed on the original velocity data [52]. For real-valued functions, *f*(*t*), *t* ∈ (−∞, ∞), the Hilbert transform is defined as a convolution of *f*(*t*) with 1/*t*:

$$H[f(t)] = \frac{1}{\pi} \int\_{-\infty}^{+\infty} \frac{f(\tau)}{t - \tau} d\tau \tag{30}$$

$$H[f(t)] = \left| H[f(t)] \right| e^{j\phi(\omega)} \begin{cases} -j \text{ } \omega > 0 \\ +j \text{ } \omega < 0 \end{cases} \tag{31}$$

where:

$$\phi(\omega) = \begin{cases} -\frac{\pi}{2} & \omega > 0 \\ +\frac{\pi}{2} & \omega < 0 \end{cases} \tag{32}$$

By performing a Hilbert transform on a real number function, the imaginary part of the function can be obtained, and the complex function can then be constructed.

Although the Hilbert transform can convert a real function into a complex function, its conditions are highly demanding, requiring the input signal to be linearly stable. In practical applications, it is difficult for most signals to meet such conditions; thus, the Hilbert transform cannot be directly used. Huang's empirical mode decomposition (EMD) can be used to decompose the signal into several

intrinsic mode functions (IMFs) and residual signals. When the residual signal is a monotonic function or lacks a maximum or minimum, the EMD ends. Here, we have:

$$\xi(t) = \sum\_{i=1}^{K} IMF\_i(t) + r\_K(t) \tag{33}$$

where ξ(*t*) is the original signal, *IMFi*(*t*) is the eigenmode function, and *rK* is the residual. Now all IMFs meet the linear steady-state requirements and can be Hilbert transforms. Finally, the transformed IMF is superimposed with the residual, that is, the original signal is converted into a complex signal. The Hilbert Huang transform is applied to the velocity field data at each timepoint as the initial data of the DMD method, and then, the standard DMD is performed.

Generally, to select a suitable mode for analysis, the modes must be sorted according to the energy. The energy of a dynamic mode can be denoted by the norm, which is defined as follows:

$$\|\|\boldsymbol{\varphi}\_{l}\|\| = \sqrt{\sum\_{i=1}^{n} \left|\phi\_{l}\right|^{2}}\tag{34}$$

where φ*i* is the velocity matrix of the *i*-th mode.

Figure 11 shows an energy spectrum of the dynamic modes. Mode 1 corresponding to a frequency of zero contains most of the energy of the flow field. Mode 2 contains the second greatest amount of energy in the flow field, and the corresponding frequency is close to the cavitation shedding frequency (35.75 Hz and 17.35 Hz) obtained in the previous section. The frequencies corresponding to mode 3 and mode 4 are close to two and three times of mode 2. Due to the strong instability of the cavitation flow field, the frequencies for mode 3 and 4 are not strictly two or three times the dominant frequency, which is acceptable. As shown in Figure 10a,b, for the near free surface case, the modes except for mode 1 contain significantly less energy than the no free surface case. This result indicates that the free surface inhibits the cavitation flow, which may explain why the length decreases for the near free surface cavitation case.

**Figure 11.** Comparison of the energy spectrum of the dynamic modes with free surface and without free surface. (**a**) Near free surface (*h* = 0.25*C*); (**b**) No free surface.

Figure 12a,b show the distributions of the real part of the dynamic modes on the mid-section plane. Mode 1 corresponds to an average flow field with a frequency of zero. The real part is negative only near the hydrofoil. The region of negative values is similar to the sheet cavitation region, and most other areas exhibit positive values. The positive and negative values alternately appear in couples for modes 2–4, reflecting the coherent structure of the turbulence. This phenomenon is due to the periodic development of the hydrofoil cavitation, which is caused as the vortex rolls up and sheds off. As the frequency increases, the scale of the coherent structure becomes smaller. Compared with the case shown in Figure 11a in the absence of a free surface, the dynamic mode scale is larger near the free surface, as shown in Figure 11b, and the coherent structure regions have a larger range. This phenomenon may arise because the turbulence generated by the hydrofoil near the free surface is not fully developed due to the influence of the free surface. Moreover, the energy contained in modes 2–4 is less when it is near the free surface.

**Figure 12.** Comparison of the dynamic modes on mid-section plane with free surface and without free surface. (**a**) Dynamic mode near the free surface (*h* = 0.25*C*); (**b**) Dynamic mode without free surface.

Reconstruction with different modes can result in flow patterns at different characteristic frequencies. In this way, the vortex motion characteristics of the flow field at different frequencies can be analyzed. Figure 13a shows the reconstructed flow field velocity vector over a single cavitation cycle near a free surface. The first column presents the numerical simulation results, and the second column displays the flow field reconstructed by all modes, which reflects the true state of the flow field. The flow field matches the numerical simulation results, demonstrating the effectiveness of the improved DMD method. Because mode 1 corresponds to a frequency of zero and it has no periodicity, it is not analyzed here. The third and fourth columns present the reconstructed flow fields for mode 2 and mode 3, respectively. It can be seen that in addition to mode 1, the coherent structure of the single-order mode is embodied as a vortex in the reconstructed flow field. Over a single period, a vortex structure is observed inside the sheet cavitation on the suction side of the hydrofoil. As the cavitation sheds off and transitions to cloud cavitation, the tails of the sheet cavitation form new vortices and move downstream. The scale of the vortex decreases as the dynamic mode increases, similar to the

mode coherent structure. When there is no free surface, as shown in Figure 13b, the reconstructed flow fields of all modes are also similar to the numerical simulation results, which express the true state of the flow field. Compared with the near free surface case, in addition to the difference from the real flow field, the reconstructed flow fields for mode 2 and mode 3 exhibit significant deviations. In the absence of a free surface, the vortex formed by the shedding cavitation appears at a position further downstream, and the scale is larger than that near the free surface. This phenomenon indicates that the free surface has an inhibitory effect on the vortex generated by the cavitation, resulting in a reduced vortex scale.

**Figure 13.** Comparison of the velocity vector of the reconstruction flow field for the near free surface case and no free surface case. (**a**) Velocity vector of the reconstruction flow field for the near free surface case (*h* = 0.25*C*); (**b**) The velocity vector of the reconstruction flow field for the no free surface case.
