1   Introduction

Cavitation is an elaborate two-phase flow, which is prevalent in hydraulic machinery. The attendant influences of cavitation, such as erosion, vibration, and noise, remarkably compromise the secure operation of hydraulic machinery (Cheng et al., 2021; Wang et al., 2021; Wu and Chen, 2016). Determining the physical mechanisms of cavitating flow is urgent. Thus, in previous decades, massive research has been conducted to study the evolution of cavitating flow concerning the growth of the sheet cavity and the shedding and collapse of the cloud cavitation (Bensow and Bark, 2010; Cheng et al., 2020; Leroux et al., 2004).

Computational fluid dynamics (CFD) is an efficient tool for gaining an in-depth knowledge of the physical systems of cavitating flow. The numerical results are normally compared with experimental data to evaluate the qualitative accuracy of the simulations (Long et al., 2018). This treatment depends on the experimental results and is unable to evaluate the simulation errors quantitatively. Verification and validation (V&V) is a common procedure for estimating the simulation error (Ayyildiz et al., 2019; Gousseau et al., 2013; ITTC, 2002; Oberkampf and Roy, 2010). The simulation error, which is defined as the difference between the simulation value and truth, can be split into numerical error and modeling error. Verification is a method for evaluating numerical errors/uncertainties, whereas validation is the assessment of modeling errors/uncertainties. Roache (1994) first proposed the grid convergence index (GCI) method for conducting V&V research. The GGI method depends on grid convergence studies with grid refinement and is the basis for the succeeding advancement of V&V methods. The later V&V research mainly concentrates on enhancing the computation of the numerical uncertainty, which entails the GCI-LN method (Logan and Nitta, 2005), the GCI-OR method (Oberkampf and Roy, 2010), and the GCI-R method (Roache, 2011). Stern et al. (2001) and Xing and Frederick (2011) further increased the application reach of the uncertainty approximation by stating the safety factor as a function. Long et al. (2017) compared the viability of these uncertainty approximation schemes in cavitating flow and learned that the approach offered by Xing and Frederick (2011) is more robust. The GCI method and its modifications are the typical procedures for V&V research, but they cannot be directly applied to the error estimation of multiscale models because the numerical and modeling errors change concurrently with grid refinement, which hinders distinguishing both error components. Large eddy simulation (LES) is a classic multiscale model that adopts spatial filtering to eliminate small-scale turbulent configuration. To determine the issue, Freitag and Klein (2006) expressed the modeling error as a power series of filter size and concurrently settled the numerical and modeling errors. Dutta and Xing (2017) and Xing (2015) further presented the H2‒5 method to assess the LES errors. To expand the range of applications, Xing (2015) also established the H2‒3 method, which regards the order of accuracy as a known quantity, such that data on three sets of grids are sufficient for computing the errors. Long et al. (2020; 2021; 2019) and Deng et al. (2021) implemented the H2‒3 method to assess the LES errors for unsteady cavitating flow simulation. The practicality of the approach has been initially confirmed. For the said error estimation approach hinged upon Richardson extrapolation, the asymptotic range condition must be met, but meeting this requirement is challenging for engineering issues. Acquiring reliable values for the order of accuracy in practical applications is difficult. Hence, scholars usually consider the order of accuracy as a known quantity to solve the problem. Some studies (Klein, 2005; Long et al., 2020; Long et al., 2019; Xing, 2015) used the order of accuracy based on macroscopic quantities to approximate the simulation errors of the local variables. The H2‒3 method is a classic instance of the said approach. However, setting the order of accuracy for the local variables in that manner is inappropriate and may lead to the miscalculation of errors. Moreover, regarding computational efficiency, researchers often used three sets of meshes in their LES V&V studies. The treatment is unable to remove the influence of the noisy data in the simulation, which lessens the integrity of error estimation outcomes. To address these problems, Eça and Hoekstra (2014) proposed a V&V approach that includes three additional error estimators, which affords a practical solution for error estimation outside the asymptotic range. Moreover, they calculated the final error results in the least-square sense to lessen the effect of the noisy data. According to the above discussion, the current paper presents a new LES V&V method and examines its viability in turbulent cavitating flow near the free surface.

Current literature examining the influence of the free surface on the cavitating flow is limited. For experimental operation, producing a stable free surface in a water tunnel is challenging. As for numerical simulations, the three-phase coupled flow creates an immense challenge to computational stability and accuracy. Regarding the above issues, Wang et al. (2016) adopted SHPB technology in experimental observations, which guaranteed the free surface static before flowing around the model. They discovered that the cavity demonstrates an asymmetric feature under the influence of the free surface. They also implemented the LES model to address the turbulent flow, and the numerical data were consistent with the experimental results. Zhou et al. (2019a; 2019b) performed numerical simulation on the same issue and reported the same phenomenon. Using SHPB technology, Xu et al. (2017b) performed experimental observations of the cavitating flow around a hydrofoil near the free surface. They learned that the cavity on the upper hydrofoil grows more gradually, and the re-entrant jet moves faster upstream. They also conducted numerical simulations to examine the effect of cavitation on vorticity distribution further. Xu et al. (2017a) studied the free surface flow around a blunt body by using numerical and experimental approaches. The effect of the free surface was similar to the results of the said studies. They later analyzed the ventilated cavitating flow around the blunt body (Wang et al., 2017). The results revealed that air entrainment promotes the creation of a larger, more stable cavity. The cavitation pattern diverges from that noted under larger submergence depth. In addition, Sun et al. (2021) examined the tip vortex cavitating flow beneath the free surface. They observed that the existence of a free surface stifles the creation and advancement of tip vortex cavitation. Studies on the cavitating flow under waves are also rare. Cheng et al. (2011) established that the shedding frequency of the cavity increases with wave height, and the presence of waves improves the re-entrant jet strength. Similar conclusions have been reached in the studies by Bal and Kinnas (2002) and Faltinsen (2005). Moreover, bubble dynamics is a valuable topic in cavitating flow, which is given increased value by the researchers. Zhang et al. (2023) developed a unified, sophisticated mathematical model to control the bubble dynamics, which considers numerous influences. Compared with the previous techniques, the theory proposed by Zhang et al. (2023) has a superior performance in replicating the bubble oscillation, migration, and collapse-induced pressure fluctuations. This theory is valuable for future studies in cavitating flow simulations (Chen et al., 2024; Wang et al., 2024; Xiong et al., 2024).

Guided by the above contents, the current paper implements LES coupled with the Schnerr–Sauer (S–S) model to examine the unsteady cavitating flow around a projectile near the free surface. A new LES V&V method is implemented to assess the simulation errors. First, the errors of the local streamwise velocity are estimated. The association between the error distribution and cavity evolution is analyzed. Second, to assess the performance of subgridscale (SGS) models, the modeling errors of two typical SGS models are compared. Third, the cavitation–vortex interaction is examined to demonstrate the influences of cavitation on the error distribution. Finally, the distribution of the monitoring points outside the asymptotic range is explained.

2   LES error estimation method

Based on previous studies (Freitag and Klein, 2006; Long et al., 2019; Xing, 2015), LES error can be split into numerical error and modeling error. According to the H2‒5 LES V&V method (Xing, 2015), the error can be assessed by the following:

where $ S_{C} $ is the numerical benchmark, and $ S_{1} $ refers to the simulation results of the corresponding case. Generally, the direct numerical simulation (DNS) results are set as the numerical benchmark. However, regarding computational cost, performing a DNS study for turbulent cavitating flow is almost impossible, so $ S_{C} $ is unknown and must be determined. $ c_{N} $ and $ c_{M} $ are undetermined coefficients. $ p_{N} $ and $ p_{M} $ are the order of accuracy for numerical and modeling errors separately. $ h^{*} $ can be derived as follows:

where $ h $ is the local grid size determined as $ h=V^{1 /3} $ ($ V $ is the volume of the local cell), $ \varDelta t $ is the time step, and $ \varDelta $ is the filter width, which is equal to the local grid size $ h $ for the implicitly filtered LES. Eq. (1) has five unknown variables and requires developing at least five equations for solving. Numerical simulations need to be performed based on at least five sets of grids. Moreover, the simulation data must satisfy the "asymptotic range", that is, the data on grids must be sufficiently fine to expand the error into a single dominant term. However, the practical application has difficulty meeting the "asymptotic range" condition. To resolve this issue, the present paper combines three additional error estimators for numerical and modeling errors to solve the problem of practical applications where the asymptotic range cannot be well satisfied. The numerical error and modeling error can be stated in the following forms:

To reduce the effect of noisy data, the current paper calculates the errors in the least-square sense based on the solutions from six sets of grids. The function to be solved is shown as follows:

where $ n $ is the number of the meshes, and subscript $ i $ is employed to identify the grids ($ i=1, 2, \cdots, 6 $).

Figure 1 presents the V&V method. First, $ \delta_{N}=\delta_{N-\mathrm{RE}} $ and $ \delta_{M}=\delta_{M-\mathrm{RE}} $ are substituted into Eq. (4). Then, $ p_{N} $ and $ p_{M} $ are determined according to the minimization of $ \Phi $. If $ p_{N} $ and $ p_{M}$ are within the asymptotic range, which is set as 0.5 ≤ p ≤ 2 according to the previous studies (Eça and Hoekstra, 2014; Phillips and Roy, 2014; Stern et al., 2006), $ \delta_{N} $ and $ \delta_{M} $ are approximated by $ \delta_{N-\text { RE }} $ and $ \delta_{M-\text { RE }} $ , respectively. If $ p_{N} $ (or $ p_{M} $) cannot meet the asymptotic range, $ \delta_{N-\mathrm{RE}} $ (or $ \delta_{M-\text { RE }} $) is dropped and replaced by $ \delta_{N-1}, \delta_{N-2} $ , and $ \delta_{N-12} $ (or $ \delta_{M-1}, \delta_{M-2} $ , and $ \delta_{M-12} $). Then, each expression is substituted into Eq. (4) to determine $ \delta_{N} $ and $ \delta_{M} $ individually. The best fit is achieved through the minimum of the standard deviation of the fit. The standard deviation is expressed as follows:

Figure  1  Diagram of V&V procedure

Download: Full-Size Img

where m is the number of unknowns.

3   Numerical methods and simulation setup

3.1   LES equations

The present paper implements the LES method with the homogeneous assumption to address the turbulent flow. The continuity and momentum equations for the vapor/water/air three-phase flow are shown as follows:

where $ u_{i} $ denotes the velocity components; $ p $ is the mixture pressure; $ \rho_{m} $ and $ \mu_{m} $ are the mixture density and laminar viscosity, respectively. They are calculated through the following equations:

where $ \alpha $ is the volume fraction of each phase; the subscripts $ l, v $, and $ \alpha $ denote the liquid, vapor, and air phases separately. In addition, the vapor and air volume fractions are determined using the following transport equations:

In Eq. (11), $ \dot{m}^{+} $ and $ \dot{m}^{-} $ are the mass transfer rates of evaporation and condensation, respectively. The current paper utilizes the S–S model to simulate the mass transfer between the water and vapor phase. The cavitation model is described below. Eq. (12) is used to track the air volume fraction throughout the domain.

The overbars in Eqs. (6) and (7) represent the filtered quantities, which are derived using the following equations:

where D is the computational domain, and G is the filter function. Moreover, the SGS stress τ_ij is presented in Eq. (14), and this term is modeled by Eqs. (15) and (16) based on the Boussinesq hypothesis.

where $ \overline{S_{i j}} $ is the rate-of-strain tensor, $ \tau_{k k} $ is the isotropic part of the SGS stress, and $ \mu_{t} $ is the SGS turbulent viscosity. The current paper utilizes the Smagorinsky – Lilly (SL) model and the Dynamic Smagorinsky–Lilly (DSL) model to compute $ \mu_{t} $. The SL model is expressed as the following equations (Smagorinsky, 1963):

The Smagorinsky constant C_s is set as 0.1. Δ is the filter width, κ is the von Kármán constant, and d is the distance to the closest wall. Because C_s is not a universal constant, the SL model may introduce numerical dissipation, which may cause the damping of large-scale pulsations. Thus, a DSL model is used to compute the constant C_s based on the resolved scales. Germano et al. (1991) and Lilly (1992) presents more detailed information about the DSL model.

3.2   Physical cavitation model

The current paper employs the S–S model (Schnerr and Sauer, 2001) to reproduce the cavitating flow, which is acquired by simplifying the Rayleigh–Plesset equation. The viability of this model has been widely proven in latest years (Hu et al., 2021; Sun et al., 2022; Zhi et al., 2022). The mass transfer rates between the water and vapor phase are expressed as follows:

where $ p_{v} $ is the vapor pressure that is set as 2 340 Pa here. $ R_{b} $ is the bubble radius determined by the following equation:

where $ N_{b} $ is the number of bubbles per unit volume. Eq. (22) is employed to link the volume fraction and bubble radius for the equation closure. $ N_{b}$ is set as a constant (10¹³) in the present paper.

3.3   Simulation setup and mesh information

Figure 2 displays the computational domain and boundary conditions. The computational domain extends 1, 0.4, and 0.6 m along the x, y, and z axes, respectively. The current paper considers half of the projectile model, which is 246 mm long. Its conical angle is 90°, and the diameter d is 37 mm. The projectile works at a depth of 17 mm underwater, and the inflow velocity U_∞ is 17.8 m/s. The outlet pressure is determined by the cavitation number and calculated as follows:

Figure  2  Computational domain and boundary conditions

Download: Full-Size Img

where p_∞ is the pressure in the open air, and a fixed pressure boundary is applied to the outlet patch; p_v is the vapor pressure set as 2 340 Pa here.

Figure 3 shows the grids around the projectile for the six sets of meshes. The Cartesian cut-cell method is employed to produce high-quality meshes, with a constant refinement ratio r = 1.2. The time step for each case shares the same refinement ratio. Table 1 presents the information for the six meshes, including the minimum face size, growth rate, and the mesh information and simulation setup for each case. Figure 4 demonstrates the y⁺ distribution for the six meshes. The maximum y⁺ for the six meshes is approximately 1, and the y⁺ value is < 1 for the most part of the projectile surface. The height of the first layer is set to 1 × 10⁻⁵ m for the six meshes. The numerical simulation is performed through ANSYS Fluent, which is based on the finite volume method. The second-order implicit scheme is employed for the transient formulation, and the bounded central differencing scheme is used for momentum equations. To resolve the discrete equations, the present paper follows the SIMPLEC algorithm, which is a pressure–velocity coupling algorithm to answer unsteady problems. Through a connection between velocity and pressure corrections, this algorithm implements mass conservation and acquires the pressure field. Regarding the model parameters, the present paper employs the S–S model to reproduce the mass transfer between the water and vapor phases. The number of vapor bubbles per unit volume is set as 1×10¹³ in the present paper, which is generally validated in the previous literature. Moreover, for the SL model, the Smagorinsky constant C_s is 0.1, and the von Kármán constant κ is 0.41. A uniform flow field is used to initialize the transient cavitating flow simulations.

Figure  3  Grids around the projectile

Download: Full-Size Img

Figure  4  Distribution of y⁺ for Meshes 1–6

Download: Full-Size Img

4   Results and discussion

4.1   Error estimation based on the new V&V method

The experimental results refer to a previous study (Wang et al., 2016). The SHPB technology (Wei et al., 2011) is employed to launch the projectile through the stress wave. The experimental setup can hasten the projectile to 30 m/s in less than 50 μs so that the free surface can stay static before launching. Figure 5 compares the cavitation patterns between the experimental and numerical results. The free surface is colored by the height in the z direction. Figure 6 compares the cavity length between the experimental and numerical results. The simulation can accurately replicate the evolution of the cavity length. At t = 1 ms, the connected cavity develops at the head of the projectile. The cavity length on the upper is almost equal to that on the lower surface, and the amplitude of the free surface is small. Then, the cavity length proceeds to extend downstream during t = 3–5 ms. Under the influence of the free surface, the cavity length starts to demonstrate an asymmetrical attribute, and the cavity length below is slightly longer than that above. The height of the free surface increases. At t = 6–9 ms, the re-entrant jet appears and slowly travels upstream, and the attached cavity sheds and moves downstream. Figure 7 exhibits the flow vectors near the projectile. The re-entrant jet is beneath the cavity, and the flow direction is opposite to the mainstream. Owing to the free surface, the velocity of the re-entrant jet above exceeds that below; consequently, the cavity above detaches and travels downstream earlier. At t = 11 ms, the shedding cavity begins to collapse. The cavity below still has some remnants, whereas the cavity above nearly collapses entirely. The shedding cavity proceeds downstream in an asymmetric manner, and the cavity below migrates a longer distance. The numerical simulations appropriately replicate the cavity evolution, including the cavity emergence, development, shedding, and collapsing downstream. The simulation using various meshes produces similar results for the attached cavity, which corresponds to the simulation results from t = 1 ms to t = 5 ms. However, as the mesh resolution rises, the simulated cloud cavitation becomes nearer to the experimental observations. For example, at t = 9 ms, the simulation results acquired from Meshes 4, 5, and 6 demonstrate richer 3D structures. Moreover, they capture the cavitation near the shoulder region, which is lacking in the results from Meshes 1 and 2. Additionally, regarding the collapse of the shedding cavity, the results acquired from Meshes 4, 5, and 6 are more thorough; the shedding cavity moves further downstream, closer to the experimental observations. Based on the simulation data, the V&V research is performed in the following contents.

Figure  5  Comparison of cavitation patterns between the experimental and numerical results

Download: Full-Size Img

Figure  6  Comparison of the cavity length between the experimental and numerical results

Download: Full-Size Img

Figure  7  Schematic diagram of the re-entrant jet

Download: Full-Size Img

Figure 8 presents the sketch of the six monitoring planes. S1, S2, and S3 are all within the cavitation region with a distance of 0.5d between each adjacent pair. S4, S5, and S6 are outside the cavitation region, and the distance between each adjacent pair is 2d. Figure 9 exhibits the distribution of dimensionless streamwise velocity along the monitoring lines at six moments. Figure 8 shows the position of the monitoring lines. With the development of cavitation, the velocity fluctuations along L1–L3 slowly intensify, specifically in the region near the projectile where cavitation arises. At t = 5 ms, the velocity opposite the main flow direction is noted near the wall due to the appearance and development of the re-entrant jet. At t = 9–11 ms, the collapse of the shedding cavity further complicates the local flow field. The velocity values slowly stabilize when moving away from the projectile. Moreover, on monitoring lines L4–L5, the velocity distribution stays stable with slight deviation between different moments. Cavitation greatly complicates the local flow field, which is accountable for the distributions of the simulation errors. Figure 10 presents the numerical and modeling errors of the flow velocity at each moment. The error results in the present paper are based on Mesh 1. The errors are normalized by the inlet velocity U_∞. Overall, the planes within the cavitation region display larger error values than those beyond the cavitation region. The extreme values are dispersed within the cavity region and the air–water interface. Moreover, as cavitation progresses, the spatial distribution of the errors extends in the spanwise and streamwise directions, and the error magnitudes rise. Cavitation complicates the local flow field, which possibly raises the simulation difficulty and subsequently intensifies the error values.

Figure  8  Monitoring plane

Download: Full-Size Img

Figure  9  Streamwise velocity distribution along the monitoring lines

Download: Full-Size Img

Figure  10  Spatial distribution of the errors on the monitoring planes

Download: Full-Size Img

At t = 1 ms, the attached cavity displays, and only the S1 plane presents a noticeable distribution of errors within the cavity region. The error values for the other planes approach zero. During t = 3–5 ms, the cavity length proceeds to increase, and the cavitating flow starts to interrupt the flow field downstream. Consequently, S2 and S3 planes also display substantial error distributions. In addition, the error magnitudes below the projectile are higher than that above the projectile in the S1 plane at 1 ms and the S2 plane at 5 ms. This phenomenon is most likely associated with the asymmetric characteristic of the cavity, as shown in Figure 11. At 1 and 5 ms, owing to the faster growth of the cavity below, the cavity below is closer to the respective plane, causing the local flow field to be more vulnerable to the cavitating flow and present greater error magnitudes. During 6–9 ms, the re-entrant jet arises and travels upstream to the head of the projectile. The intense shear flow further upsets the local flow field, leading to a wider range of error distribution and rising error magnitudes. Lastly, the shedding cavity asymmetrically collapses at 11 ms. The error distribution also demonstrates an asymmetrical feature at the S1, S2, and S3 planes, where the error distribution below the projectile displays a greater area and magnitude. This finding may be ascribed to the intense velocity variations caused by the collapse of the nether cavity at that moment, causing greater local errors. The cavity above has nearly collapsed completely, and the local flow field stabilizes compared with that below the projectile, resulting in slighter error magnitudes. The cavity evolution has a substantial effect on the error distribution. Moreover, the present V&V approach is strong within the cavity and free surface region. The results verify the applicability of the V&V method in cavitating flow simulation.

Figure  11  Position of the cavity and planes

Download: Full-Size Img

4.2   Modeling error comparison between the SL and DSL models

This section further examines the applicability of the V&V method for unsteady cavitating flow. Figure 12 compares the modeling errors between the SL and DSL models. The DSL model shows a distribution of modeling errors similar to that of the SL model. Greater error magnitudes for the two models are concentrated around the cavity region and the water–air interface. However, the key difference is that the error magnitudes of the DSL model are smaller than those of the SL model. Figure 12 shows that during t = 1–3 ms, the error magnitudes of the DSL model are evidently less than that of the SL model within the cavity region. During t = 6–11 ms, the error results show that the DSL model can more correctly address the turbulent variations generated by the re-entrant jet and collapse of the shedding cavity. Especially at 11 ms, the DSL model demonstrates a remarkably smaller error distribution compared with the SL model. The DSL model attains a higher accuracy in reproducing the unsteady cavitating flow. This finding can be ascribed to the numerical dissipation induced by the utilization of the constant coefficient C_s, as shown in Eq. (18) in the SL model. However, the DSL model computes the coefficient based on the resolved scale, which suggests that the DSL model dynamically alters the coefficient based on the flow field, making it more appropriate for the local flow conditions. This conclusion is commonly verified in previous studies (Cheng et al., 2020; Ekman et al., 2021; Ghasempour et al., 2014), and the error results here also demonstrate the point. The present error estimation procedure can capture the reduction in modeling error owing to model enhancements, which demonstrates its effectiveness for unsteady cavitating flow concerns.

Figure  12  Comparison of the modeling error between the SL and DSL models

Download: Full-Size Img

To compare the model's capability to resolve the turbulent fluctuations, Figure 13 presents the resolved Reynolds stress between the SL and DSL models. $ \overline{u^{\prime} u^{\prime}} $ is the Reynolds normal stress. τ denotes the Reynolds shear stress, which is determined as $ \tau=\sqrt{\left(u^{\prime} v^{\prime}\right)^{2}+\left(u^{\prime} w^{\prime}\right)^{2}+\left(v^{\prime} w^{\prime}\right)^{2}} $. $ \overline{u^{\prime} u^{\prime}} $ and τ are all normalized by $ U_{\infty}^{2} $. The Reynolds normal stress exhibits a definite distribution near the projectile wall, whereas the large Reynolds shear stress is concentrated within the shedding cavity region. The velocity fluctuations caused by the re-entrant jet are most likely the source of the Reynolds normal stress. The shedding cavity collapse is the major factor causing an increase in the Reynolds shear stress. The DSL model can resolve a broader range of Reynolds stress compared with the SL model. The DSL model can better capture the unstable variations due to the re-entrant jet and the cavity collapse, which is also confirmed by the error results presented in Figure 12. In addition, the more turbulent fluctuations can be captured with the refinement of the grids. Mesh 5 has more elements than Mesh 4, but few enhancements are noted in Figure 13. Moreover, Figure 5 reveals that Mesh 4 can appropriately model the cavity evolution, specifically the cloud cavitation and collapse of the shedding cavity. The setup for Mesh 4 is adequate for performing the LES simulation in the present paper. The following discussion is based on Mesh 4 considering the simulation accuracy and computational consumption.

Figure  13  Comparison of the resolved Reynolds stress between the SL and DSL models

Download: Full-Size Img

4.3   Data outside the asymptotic range

Satisfying the asymptotic range condition is challenging for the cavitating flow simulations. To prevent that from causing the abandonment of the error estimation, the current paper suggests a new V&V method. The above contents confirm its reliability for the cavitating flow simulation. This section primarily explains the simulation data beyond the asymptotic range. More than 10⁴ monitoring points are marked on each monitoring plane, as shown in Figure 8. Then, the streamwise velocity at each monitoring point is obtained to perform the V&V research. The V&V results are presented in the above sections. Figure 14 shows the simulation data beyond the asymptotic range at each moment. The black dots denote the monitoring points where the data meet the asymptotic range condition, whereas the red dots represent the monitoring points where the data are outside the asymptotic range. The monitoring points where the data are beyond the asymptotic range mainly disperse around the cavity region. Some such points are also distributed near the free surface, but the number is substantially smaller. As the cavity evolves, the number of these points rises, and their distribution area slowly increases, specifically during the re-entrant development phase and the shedding cavity collapse phase. Nonetheless, on the S4, S5, and S6 planes, which are beyond the cavity region, a sparse distribution of these points is consistently observed. Thus, the local complex flow introduces additional challenges in correctly determining the order of accuracy, causing the inability to satisfy the asymptotic range condition.

Figure  14  Monitoring points where the data are beyond the asymptotic range (The red dots represent the monitoring points where the data are outside the asymptotic range)

Download: Full-Size Img

The order of accuracy $ p_{N} $ and $ p_{M} $ are of specific interest in the existing V&V studies and also the concern in the current paper. Thus, the contours of the estimated $ p_{N} $ and $ p_{M} $ at typical X positions are added. Figure 15 presents the spatial distribution of $ p_{N} $ and $ p_{M} $ on the monitoring planes. The magnitudes of $ p_{N} $ and $ p_{M} $ are small near the free surface, resulting in greater error magnitudes in that region. This occurrence becomes more evident with the advancement of the cavitating flow, which is in line with the analysis above. In the cavitation region, the distribution of $ p_{N} $ and $ p_{M} $ is complicated and exhibits apparent pulsations. As the cavitating flow proceeds, these pulsations intensify, and the range of fluctuations grows, leading to remarkable oscillations in the local error values. By contrast, on the S4–S6 planes, the distribution of $ p_{N} $ and $ p_{M} $ is more constant, which facilitates the lower and steadier error results on the planes.

Figure  15  Spatial distribution of $ p_{N} $ and $ p_{M} $ on the monitoring planes

Download: Full-Size Img

Figure 16 presents the number of monitoring points where the data are beyond the asymptotic range at the typical moment. The horizontal axis represents the various planes in Figure 8, whereas the vertical axis n denotes the number of monitoring points where the data are beyond the asymptotic range. As time advances, the value of n substantially rises on the S1, S2, and S3 planes, but its increase on the other planes is gentler, almost unchanged. The increasingly strong cavitating flow hinders accurately determining the order of accuracy, which most likely leads to the sustained growth of n. Moreover, the n values of the DSL model are evidently less than those of the SL model within S1–S3 planes. Figure 17 also presents the evidence for the viewpoint. Figures 17(a) and 17(b) exhibit the n values for each plane during t = 1–11 ms, and Figure 17(c) demonstrates the total number of monitoring points where the data are beyond the asymptotic range for the six planes. Particularly, during the shedding cavity collapse phase, the n value of the DSL model is considerably less than that of the SL model. This result may be ascribed to the superior capability of the DSL model in addressing unstable fluctuations. Therefore, the simulation data based on the DSL model are adequately acceptable for approximating the order of accuracy, influencing the reduction of the monitoring points where the simulation data are outside the asymptotic range.

Figure  16  Number of the monitoring points where the data are beyond the asymptotic range

Download: Full-Size Img

Figure  17  Comparison of the number of monitoring points beyond the asymptotic range between the SL and DSL models

Download: Full-Size Img

The cavitating flow prevents the simulation data from satisfying the asymptotic range, which limits the implementation of V&V in cavitating flow simulations. The present paper presents a new V&V method for LES, which supplies a viable solution for the error approximation of simulation data that cannot satisfy the asymptotic range. Sections 4.1 and 4.2 discuss that the new V&V method performs robustly within the cavitation region, where the simulation data are beyond the asymptotic range. Moreover, the V&V method can capture the decrease in modeling error owing to model enhancements. This paper offers a solution for the expansion of V&V application in cavitating flow simulations.

4.4   Cavitation–vortex interaction and V&V results

To investigate the correlation between cavitation and error distribution further, the influences of cavitation on the flow structure are explained. Figure 18 exhibits the vortex structure identified by the Q-criterion (Q = 10 000 s⁻²) around the projectile. The Q-criterion is expressed as follows:

Figure  18  Vortex structure around the projectile

Download: Full-Size Img

where Ω is the vorticity rate, and S is the strain rate. The vortex structure is primarily distributed within the cavity region. During t = 1–3 ms, the vortex structure positioned at the head of the projectile steadily develops downstream with the increase of the cavity length. The vortex structure arises quickly within the cavity region, complicating the local flow field and enhancing the simulation errors. Beyond the cavity region, fewer vortex structures are dispersed, and the stable flow field may lead to smaller local simulation errors. After 6 ms, the re-entrant jet emerges and travels upstream, prompting powerful vortex motion in the local region. The vortex structure begins to demonstrate an elaborate 3D characteristic and extends downstream with the shedding cavity migration. At 11 ms, the collapse of the shedding cavity further initiates the creation of the largescale vortex structures downstream. They establish complex dynamics in the surrounding flow field, leading to fluctuations in velocity and flow patterns. This situation most likely increases the simulation difficulty and the local error magnitudes.

Figure 19 exhibits the vorticity distribution on the monitoring planes. The vorticity is concentrated within the cavity region. The vorticity at the same moment slowly decreases on each plane when traveling downstream, owing to the decreased influence from the cavitating flow. Moreover, due to the asymmetric characteristic of the cavity, the vorticity develops first below the projectile on the corresponding plane, such as S1, S2, and S3 planes at 1, 5, and 9 ms, respectively. This outcome results in pronounced fluctuations in the local velocity field, which may be associated with the asymmetric distribution of errors on the corresponding plane. Moreover, Figure 20 shows the influences of cavitation on the flow field through the vorticity transport equation defined as follows:

Figure  19  Vorticity distribution on the S1, S2 and S3 plane

Download: Full-Size Img

Figure  20  Contours of the vortex stretching, dilatation, and baroclinic torque terms

Download: Full-Size Img

The stretching term $ (\boldsymbol{\omega} \cdot \nabla) \boldsymbol{V} $ denotes the effects of the stretching and tilting of vortices. The dilation term $ \boldsymbol{\omega}(\nabla \cdot \boldsymbol{V}) $ represents the influence of fluid compressibility on the vorticity distribution. The baroclinic torque term $ \frac{\nabla \rho_{m} \times \nabla p}{\rho_{m}^{2}} $ is caused by the misalignment between the pressure and density gradients. The last term considers the effects of viscous diffusion, which is much slighter than the three other terms (Ji et al., 2017; Ji et al., 2015). The last term is ignored in the subsequent discussion. Figure 20 shows that the stretching term dominates among the three terms, which presents an evident distribution near the cavity region and water– air interface. Specifically at 9 ms, the re-entrant jet triggers an abrupt change of the local velocity gradient, amplifying the magnitudes of the stretching term. At 11 ms, the collapse of the shedding cavity further intensifies the local velocity fluctuations, which substantially improves the stretching term. The dilation term principally exists within the cavity area and has nearly no distribution near the water–air interface. Its distribution area increases as the attached cavity develops. The baroclinic torque and stretching term have similar distributions, and the former has smaller magnitudes. In addition, the stretching and baroclinic torque terms have positive values at the water–air interface, which means these two terms boost the generation of the local vorticity near the free surface. Inside the cavity, the signs of the three terms alternate regularly, demonstrating the vorticity within the cavity region changes dramatically.

5   Conclusions

The current paper presents a new LES V&V method and validates its applicability for the unsteady cavitating flow near the free surface. To perform the V&V research, LES coupled with the S-S cavitation model is employed in simulating the cavitating flow around the projectile near the free surface. The main conclusions are as follows:

1) The cavitating flow hinders the simulation data from meeting the asymptotic range, which limits the application of V&V in cavitating flow simulations. The new V&V method offers a viable solution for the error estimation of simulation data that cannot satisfy the asymptotic range. The V&V results reveal that larger error magnitudes are mainly distributed around the cavity region and free surface. The error values outside the cavitation region approach zero. In addition, as the cavity evolves, the increasingly intense cavitating flow poses an increasing challenge to the simulation accuracy, most likely leading to larger error magnitudes and a wider distribution area. Moreover, the error distribution occasionally demonstrates an asymmetric characteristic, which may be associated with the asymmetric morphology of the cavity at the corresponding moments.

2) The DSL model demonstrates a distribution of modeling errors similar to that of the SL model, but the modeling error magnitudes of the former are evidently smaller than those of the latter. The DSL model realizes higher accuracy in simulating the unsteady cavitating flow compared with the SL model. That viewpoint has been extensively verified in previous studies. The present V&V method can capture the reduction in the modeling errors owing to model enhancements, further exhibiting its applicability in cavitating flow simulations. Moreover, the DSL model can resolve a wider range of Reynolds stress compared with the SL model, specifically within the regions influenced by the re-entrant jet and shedding cavity, which validates the advantages of the DSL model in solving turbulent fluctuations.

3) The monitoring points where the simulation data are beyond the asymptotic range are mainly distributed near the cavity region, and the number of such points increases as the cavitating flow becomes increasingly intense. The local complex flow introduces additional challenges in accurately determining the order of accuracy, most likely causing the inability to meet the asymptotic range condition. Moreover, the DSL model has fewer such points than the SL model, which may be ascribed to the superior capacity of the DSL model in simulating turbulent fluctuations. The simulation data based on the DSL model is adequately fair for estimating the order of accuracy, resulting in the decrease of such points. The model enhancements can strengthen the ability to simulate the local complex flow, hence easing the "asymptotic range" issue.

4) The unsteady cavitating flow around the projectile near the free surface is accurately reproduced, involving the attached cavity appearance, the occurrence of a reentrant jet, and the collapse of the shedding cavity. Moreover, cavitation has a substantial influence on the flow structure. The vortex structure is concentrated in the cavity region. The re-entrant jet and shedding cavity collapse are the principal sources of the vorticity motions. The vortex motions slowly weaken away from the cavity region. The vorticity transport equation shows that the stretching term and baroclinic torque term have similar distributions around the cavity and free surface region and these two terms promote the creation of the local vorticity near the free surface. The dilation term principally exists within the cavity and has almost no distribution near the water–air interface.