We study the cavitating flow over a backward facing step with an incompressible polydisperse cavitation model. The model can predict experimental observations for this flow reasonably well, including the shedding cloud characterized by the condensation front, cavity length, void fraction, and shedding frequency. All model variations produced shedding cavities, but the turbulence model and grid resolution are essential for better predictions, with delayed detached eddy simulation (DDES) performing better than Reynoldsaveraged Navier–Stokes. Quantities, such as pressures at key points, maximum void fraction location, and shedding frequency, are mildly sensitive to those factors. Finer DDES grid resolution, crucial to resolve small vortices where cavitation occurs in their low pressure cores, improves predictions. Since a fully incompressible model produces a condensation front that follows well the experimental trends, it is concluded that compressibility is not a necessary condition for the formation of a condensation front. Consequently, the speed of sound in the mixture does not appear to play an important role in the front formation and evolution. The polydisperse nature of the model allows prediction of the bubble size distribution. Small bubbles concentrate on the downstream section of the cavity, where cavity collapse is strongest and bubble fission is most intense, while larger bubbles reside near the step where the flow is milder. The condensation front is a moving source of vorticity for the liquid phase where the “compressibility,” in the sense of mixture density changes due to void fraction changes, and baroclinic effects are significant, but the buoyancy effect is negligible.
I. INTRODUCTION
Cavitation occurs where the pressure drops below the vapor pressure. Due to its detrimental effects to machinery operation, it has attracted significant interest from researchers to develop reliable predictive models. The widely used homogenous mixture model has been successfully applied to several cavitation flows.^{1–4} Hybrid Eulerian–Lagrangian models have been in active development over several years.^{5–9} Recently, Li and Carrica^{10} developed a cavitation model based on an Eulerian approach where bubbles are described by a bubble number density distribution function. Though the model is on its infant stage, it has the potential to expand over other existing approaches by considering and predicting the bubble size distribution.
As with any model aspiring to be useful for engineering applications, iterative validation and improvements are mandatory to gain robustness and confidence in the results. The experimental data available as validation cases for cavitation models typically lack twophase information, such as bubble size distribution and void fraction, and therefore miss valuable insight into guide modeling of important multiphase flow phenomena like evaporation/condensation and bubble breakup and coalescence. As more work has been done to measure the bubble nuclei size distribution, void fraction in cavities, noncondensable bubbles in the wake, and other twophase flow quantities,^{11–17} the database for the validation of cavitation models has considerably expanded to the benefit of modelers. Recently, Bhatt et al.^{18} measured the cavitating flow behind the backward facing step (BFS) using xray densitometry and provided abundant data for the void fraction of the unsteady flow. The BFS flow is relevant since it features a turbulent shear layer, incoming turbulent boundary layer, and flow recirculation and reattachment, which are very challenging to predict for the computational models even without cavitation.
In the flow over a BFS cavitation, inception occurs first in the shear layer as streamwise “braids” (vortex structures) are stretched by larger spanwise vortices.^{19} The experiments done by Agarwal et al.^{20} on a similar BFS as in Bhatt et al.^{18} clearly showed 1–2 mm in diameter and 5–7 mm in length cavitating structures for all conditions. Discrete elongated bubbles are also detected, providing an indication for the inception region. As the cavitation number is decreased, the bubbles start to fill the region of the shear layer. Further reducing the cavitation number results in partial cavitation where cavitating structures form behind the step that behave as unsteady shedding clouds. The selfsustained oscillation of the developed partial cavity shows some similarities but of different nature to the transient cavity shedding observed in the twisted foil,^{21} which is caused by a reentrant jet. Historically, two main mechanisms have been proposed for the cavity shedding: reentrant jets and condensation shocks.^{11,22,23} Recently, Zhang et al.^{24} proposed the collapseinduced pressure wave caused by the shed cloud as another mechanism for cavity shedding based on the experimental observations. Such pressure wave can disturb the growing cavity and travels faster than the condensation shock in the mixture flow. In addition, it does not result in an apparent discontinuity in void fraction. Bhatt et al.^{18} proposed that the shedding (and collapse) of partial cavitation behind the BFS is driven by the condensation shocks. At lower cavitation numbers but before supercavitation, a cavitating flow regime is identified by Bhatt et al.,^{18} where cavities and bubbles fill a large region behind the step and the flow is quite statistically steady.
Though in terms of the jump condition, the moving front closing the cavity is analogous to the propagating shock waves in a bubbly mixture as reported by Ganesh et al.,^{11} we prefer the term condensation front instead of condensation shock. This is because shocks are most often related to fluid compressibility and the speed of sound of the liquid/vapor mixture; a compressible solver is necessary to capture these shocks. In this paper, we use mostly the term condensation front when referring to the phenomena of rapid condensation observed in the flow over the backward facing step, to prevent confusion with traditional shocks in compressible flow. Exceptions are made when discussing the speed of the front obtained from analytical shock relationships and the Mach number in the mixture, to be consistent with the terminology used in the literature.^{18,25}
Condensation fronts in cavitating flows have been observed numerically and experimentally, expressed as a moving front associated with significant condensation phase change.^{11,22,25,26} The phenomenon is quite different from shocks produced by the collapse of cavities, reported as a highly compressible process involving high pressure differences. Note that those collapseinduced shocks can trigger the condensation front formation,^{25,27} but they are not a necessary condition for the formation of condensation fronts.^{24,25} In the case of condensation fronts, the pressure differences across the jump are very low (within a few kPa), as reported in the literature.^{11,25} However, this condensation front occurs with high void fraction in the flow and alters the subsequent dynamics of the shedding process for the partial cavity. Due to the high void fraction involved, the speed of sound decreases significantly.
The vast majority of numerical simulations presenting condensation fronts in the literature employ compressible approaches.^{25–27,30} In the limited literature using incompressible solvers, CoutierDelgosha^{31} showed the capture of a condensation front. Though not shown here, the authors also performed simulations of the 2D BFS discussed in this paper using the commercial software FLUENT 22R2 with the incompressible mixture flow solver and the Zwart–Gerber–Belamri cavitation model.^{32} The condensation front was captured, though shedding frequency and void fraction exhibit large differences with the experimental data. The incompressible twofluid cavitation model used in this work is also capable of producing a condensation front, indicating that its formation is in fact not related to a process dominated by traditional shock waves associated with large pressure jumps (such as collapseinduced shocks), as discussed earlier.
In this paper, we employ the polydisperse cavitation model developed by Li and Carrica^{10} to study the cavitating flow over a BFS focusing on both model validation and flow analysis. A short description of the cavitation model is presented first, and then simulations under different conditions are introduced, along with detailed comparison with the experimental data, flow analysis, and discussion.
II. MATHEMATICAL MODELING
A. Mass transfer model
B. Bubble coalescence model
C. Bubble breakup model
D. Nuclei treatment
It is known that cavitation can occur from free nuclei in the flow (homogeneous nucleation) or surface nucleation sites due to imperfections on solid surfaces (heterogeneous nucleation). The cavitation model developed by Li and Carrica^{10} naturally incorporate these two kinds of nucleation. As a result, the nuclei size distribution and nucleation site density are input parameters for the model and obtained from experiments or estimations.
E. Limitations of the cavitation model
As with all models, the model presented herein has the following limitations resulting from necessary assumptions and simplifications:

The model cannot handle 100% void fraction (pure vapor), as the twofluid model in which is based needs a continuum liquid.

The model cannot capture pressure waves in the flow since bubbles and liquid are assumed incompressible.

The model is not suited for the prediction of cavitation inception, which is sensitive to bubble dynamics. From Eq. (11), the model would predict evaporation every time the pressure drops below vapor pressure, regardless of the duration of the low pressure event.

The accuracy of the bubble/liquid interaction models decreases as the void fraction increases.

The model is not able to predict the nuclei size distribution in noncavitating free stream regions if breakup and coalescence occur there.
Further model development could resolve some of these limitations.
III. NUMERICAL METHODS
The cavitation model was implemented in the general purpose incompressible computational fluid dynamics (CFD) code REX. REX employs multiblock structured curvilinear bodyfitted grids with dynamic overset capabilities.^{39} A singlephase level set approach is used for free surface modeling.^{40} Different RANS turbulence models with their relevant detached eddy simulation (DES) versions and boundary layer transition models were implemented for turbulence modeling.^{41} In addition, various large eddy simulation (LES) models with algebraic wall models are available for high fidelity simulations. The governing equations are discretized using a control volume/finite differences approach on a collocated grid in space. The convection terms are discretized with higherorder linear upwind biased or highorder nonlinear schemes including flux/slope limiter like total variation diminishing (TVD) schemes (up to fourth order) and essentially nonoscillatory/weighted essentially nonoscillatory (ENO/WENO) schemes. The diffusion terms use secondorder central differences. A projection method is utilized for pressure–velocity coupling including disperse phase effects,^{35,42} leading to a Poisson equation for the pressure and a linear sparse system of equations that is preconditioned and solved iteratively with Krylov solvers available in the PETSc package.^{43}
The full nonlinear system of governing equations is solved using a partitioned approach in which each variable is solved sequentially within a fixedpoint iteration for each time marching step. For the cavitation terms, time splitting is employed to enable large time steps with strong cavitation sources, though the temporal order is reduced to first order. The scheme is shown to be stable with the dynamic time step technique used in this work.
IV. SIMULATION DETAILS
Bhatt et al.^{18} studied the cavitating flow over a BFS. The experiments were performed in a recirculating cavitation channel with equal height and span of 78.7 mm. The step height $H$ was 12.1 mm. Data include measurements of the void fraction using xray densitometry under different Reynolds and cavitation numbers. Three types of cavitation regimes were observed: developing, selfsustained shedding, and close to supercavitation. As mentioned by the authors, the developing cavitation regime exhibits very low void fraction and large uncertainty exists for the data. The selfsustained shedding regime is chosen for the simulation, of the most interest due to its transient nature and best suited for the cavitation model developed by Li and Carrica.^{10}
Figure 1 shows the computational domain and grid with refinement in the region of interest. The total number of grid points is 3.5 M for the coarse grid. We perform a grid study with the refinement ratio of $ 2$ in each direction to investigate grid independence of the cavitation model, with the fine grid containing 28.5 M grid points. The incoming flow into the step is a fully developed turbulent boundary layer. In this study, we neglect explicit pressure and velocity fluctuations in the incoming turbulent boundary layer, assuming that the effects on the naturally highly turbulent cavitating flow on the BFS are small. In order to obtain the proper inlet mean velocity profile, the scaled data from a similar experiment by Agarwal et al.^{44} are employed and the velocity is imposed as a fitted polynomial. Periodic boundary conditions are applied in the spanwise direction. The inlet mean velocity for the simulation is 10 m/s, resulting in a Reynolds number of $ 1.21 \xd7 10 5$ based on the step height.
The experiments do not report the nuclei size distribution; thus, the generic size distribution suggested by the authors^{10} is employed. Note that heterogeneous surface nucleation has negligible effect on this cavitating flow. To study the effects of different modeling terms in the twophase model, we perform a series simulation as listed in Table I. Most of the cases uses the coarse grid to observe trends while minimizing the computational cost. Overall, there are three categories of studies: grid independence, turbulence model effects, and sensitivity to the cavitation number. Other details and settings for the cavitation model are the same as in Li and Carrica,^{10} unless explicitly stated in the table. 15 bubble size groups with radii logarithmically distributed between [0.01 and 2 mm] are used in all the simulations.
Case name .  Conditions . 

C_R  RANS with coarse grid 
C_D  DDES with coarse grid 
M_D  DDES with medium grid 
F_D  DDES with fine grid 
$ C _ \sigma L$  C_D case with $\sigma $ decreased 0.04 (−2 kPa reference pressure) 
$ C _ \sigma H$  C_D case with $\sigma $ increased 0.04 (+2 kpa reference pressure) 
Case name .  Conditions . 

C_R  RANS with coarse grid 
C_D  DDES with coarse grid 
M_D  DDES with medium grid 
F_D  DDES with fine grid 
$ C _ \sigma L$  C_D case with $\sigma $ decreased 0.04 (−2 kPa reference pressure) 
$ C _ \sigma H$  C_D case with $\sigma $ increased 0.04 (+2 kpa reference pressure) 
The experiments define the cavitation number based on the pressure measured on a location upstream of the step beyond the simulation domain. However, the pressure was also measured downstream of the step where there is no cavitation (see point 4 in Bhatt et al.^{18}) The pressure difference between these two locations is around $ 5.25 \xd7 10 3 \u2009 Pa$ for the case studied in this work, so the cavitation number based on the pressure at the location downstream of the step and the incoming velocity is $\sigma =0.415$. This corresponds to $\sigma =0.72$ in the experiment, with a reported standard deviation of 0.03–0.04. In the CFD simulations, the twophase flow solver forces the total pressure to satisfy the cavitation number observed in the experiment at every time step.
V. RESULTS
An animation from the simulations using the fine DDES grid is shown in video I in the supplementary material. For analysis, solutions are averaged in time and in the spanwise direction for all cases. In addition, the time window for the averages covers seven or more shedding periods except for the case with low cavitation number.
A. Flow field
Figure 2 presents a part of the time history of the total volume of vapor in the computational domain defined as $ V = \u222b \Omega \alpha d d v$, where $\Omega $ indicates the whole domain. After an initial transient, the vapor volume evolution is quasiperiodic exhibiting different frequencies and amplitudes for each condition, indicating that all cases predict an unsteady shedding vapor cloud. The amplitude of the total vapor volume indicates higher void fraction or spatial extent of the vapor cloud, while the period is controlled by shedding. As the grid is refined, the total vapor volume decreases slightly and no much difference is observed between medium and fine grid results. While period and amplitude are affected by grid refinement, grid convergence is observed for the mean total volume (see Table II). RANS predicts more vapor in the flow than DDES for the coarse grid. In addition, the magnitude and frequency of the vapor volume fluctuations are strongly affected by the cavitation number. For the low cavitation number case $ C _ \sigma L$, a large, long cavitation cloud forms. The shedding period is discussed later in this paper.
Case name .  $ L C$ .  $ L M$ .  $ L R$ .  $ V \xaf$ ( $ m 3$) . 

C_R  9.37  3.20  9.1  1.31 × 10^{−5} 
C_D  7.93  3.01  8.0  1.06 × 10^{−5} 
M_D  7.01  2.92  7.4  8.12 × 10^{−6} 
F_D  6.82  2.87  7.2  7.91 × 10^{−6} 
$ C _ \sigma L$  10.62  2.43  8.0  5.82 × 10^{−6} 
$ C _ \sigma H$  6.67  3.25  6.9  1.31 × 10^{−5} 
Case name .  $ L C$ .  $ L M$ .  $ L R$ .  $ V \xaf$ ( $ m 3$) . 

C_R  9.37  3.20  9.1  1.31 × 10^{−5} 
C_D  7.93  3.01  8.0  1.06 × 10^{−5} 
M_D  7.01  2.92  7.4  8.12 × 10^{−6} 
F_D  6.82  2.87  7.2  7.91 × 10^{−6} 
$ C _ \sigma L$  10.62  2.43  8.0  5.82 × 10^{−6} 
$ C _ \sigma H$  6.67  3.25  6.9  1.31 × 10^{−5} 
The mean streamwise velocities for C_D and F_D are shown in Fig. 3. A large recirculation is present, with a secondary corner vortex captured with both grids. The other cases predict similar vortices on the mean flow, but with different reattachment lengths, as shown in Table II. The experimental reattachment length for a singlephase flow is approximately 5.4H, while the reattachment length in cavitating flow was not measured. The longest reattachment is 9.1H predicted by C_R. As expected, as the cavitation number increases, the reattachment length decreases, approaching the value for noncavitating flow. Conversely, the reattachment length increases as the cavitation number decreases. Note that the uncertainty in the reattachment length is larger for $ C _ \sigma L$ as statistics are poorer due to the large shedding period and also due to lower grid resolution at the cavity tail, which extends into the coarse grid region.
Figure 4 presents the mean void fraction for the cases in Fig. 3. Overall, C_D predicts higher void fraction and a larger cavity, though the main shapes are similar. Cavitation initiates at the lip of the step, with much lower void fraction in that area with the finer grid. For the far wake, 7H or more downstream from the step, most of the bubbles locate away from the bottom wall, primarily due to the cavitating shedding cloud as discussed later. A similar trend is observed in the other cases. To characterize the cavity, the cavity length $ L C$ is defined as the farthest distance to the step where the void fraction reaches $ \alpha d = 0.05$ ( $ \alpha d = 0.06$ for $ C _ \sigma L$ due to poor statistics). $ L M$ is the distance of the maximum void fraction location to the step. As shown in Table II, the length $ L M$ is close to 3H for the cases with the same cavitation number, indicating the maximum void fraction location is not strongly affected by grid refinement or the turbulence model. However, the cavity length varies significantly, with RANS C_R showing the longest cavity length. Figure 5 compares the model predictions with the experimental data. The maximum void fraction location $ L M$ is predicted satisfactorily, even for different cavitation numbers. The cavity length predicted by F_D shows excellent agreement with the experimental data. Using a RANS turbulence model increases the cavity length. For the low cavitation number case, the overshoot is considerable, possibly affected by poorer statistics due to the large period between cloud shedding.
The mean void fraction profiles along the y coordinate (vertical lines) at different distances from the step are presented in Fig. 6. As expected and in line with results shown in Figs. 4 and 5, C_R predicts higher void fraction in most of the regions compared with the experiments and the C_D case. All the cases overpredict the void fraction. Li and Carrica^{10} show that the turbulence model is crucial to predict the unsteady shedding cloud in a twisted foil. In this work, that conclusion is extended, showing that the turbulence model is important for the prediction of void fraction. Higher grid resolution results in better results, though F_D still slightly overpredicts the void fraction at all stations. The result converges over different grids with respect to the void fractions at acquisition locations. The solutions for F_D are employed in the analysis of the flow.
B. Bubble size distribution
The main contribution of the cavitation model presented herein is that it can predict the bubble size distribution. The mean void fractions for different group bubble sizes are presented in Fig. 7. The maximum void fraction location moves toward the step and closer to the bottom wall for larger bubbles. Note that the void fraction magnitude is quite different between groups and that the total void fraction is dominated by the large bubbles. Though the void fraction is low at $ 10 \u2212 5 \u2013 10 \u2212 4$ for the smallest bubbles, the number density of these bubbles is very high as the bubble volume is very small. The location of the highest void fraction in group 15 is close to that for the total void fraction shown in Fig. 4, as the largest bubbles comprise most of the vapor volume. The smallest bubbles have relative low void fraction there but are present further downstream. This is because small bubbles are produced by bubble fission caused by the strong collapse of the cavity. This accumulation of bubbles of different sizes at different locations indicates that the evaporation and condensation rates can significantly vary in the flow even with the same pressure since the source is strongly dependent on the bubble size, see Eq. (11).
The mean bubble size, represented by the Sauter mean diameter $ d 32$, and the average diameter for all bubbles $ d 10$ are shown in Fig. 8. $ d 10$ exhibits a maximum of $ d 10 \u223c 120 \u200a \u2009 \mu m$ at the bottom wall near the step ( $ x = 1.5 H$), where the residence time is high and the bubbles are exposed to low pressure for longer periods. In the area where the total void fraction is very high, $ d 10$ is around 80 $ \mu m$. There is a long tail with relatively high $ d 10 \u223c 60 \u200a \mu m$, which corresponds to small vortex cavitation downstream of the main cavity (discussed later in reference to Fig. 14). In addition, very small sizes are observed on the fringe of the tail, due to bubble fission during strong condensation. Large values of $ d 32 \u223c 2 \u2009 mm$ are found very close to the step, consistent with the large bubbles residing for long periods observed in the experimental video of Bhatt et al.^{18} In the regions of highest void fraction, $ d 32$ is around 1.5 mm. The mean size becomes smaller downstream of the step due to strong condensation in those locations. In the region around $ x = 7 H$, both $ d 10$ and $ d 32$ are relatively small, corresponding to the tail of the large cloud. The larger values farther downstream are probably caused by the occasionally large cavities shed from the main cavity and transported with the flow.
The bubble size distribution was extracted at key representative monitor points shown in Fig. 9. As presented in Fig. 10, the profiles in locations close to the step (indicated by M1, M2, and M3) show two distinguishable slopes with separation at around $ 50 \u2009 \mu m$. The profiles shift upward for M1, M2, and M3 indicating an increase in the void fraction. For M4, M5, and M6, the bubble size distributions do not change significantly, showing a bit lower concentration of larger bubbles at M6. For M7, located in the wake at $ x \u223c 7 H$, the slope becomes steeper for bubbles smaller than $ 200 \u2009 \mu m$, indicating the presence of large amounts of smaller bubbles, consistent with the small $ d 10$ shown in Fig. 8. The accumulation at the largest bubble group is not as considerable as with the other three monitor points in the same panel of Fig. 10, consistent with the decrease in $ d 32$ shown in Fig. 8 for that location. All monitor points exhibit a small accumulation at $ R = 2 \u2009 mm$, suggesting a buildup of the largest bubbles.
C. Total pressure
The total pressure was investigated at three locations where there is experimental data available, as indicated in Fig. 11. The mean values for all cases are presented in Table III. All cases yield good predictions for the mean total pressure at points $ p 1$ and $ p 2$, with the maximum error at less than $ 5 %$ and $ 10 %$, respectively. Overall, the fine grid produces better results for the DDES cases, but not for all measured quantities. Larger errors are present for $ p 3$, located close to the reattachment point as shown in Fig. 11. The mean pressure contours show considerable pressure gradients in that region, and consequently, an accurate prediction of the reattachment length is important to predict the pressure in $ p 3$. All cases underpredict the mean pressure at $ p 3$ with the RANS case C_R showing the largest error. This underprediction may be due to all cases overpredicting the mean reattachment length, as shown in Table II, though the lack of experimental data prevent us from testing this assumption.
Case name .  $ p 1$ (Pa) .  Error (%) .  $ p 2$ (Pa) .  Error (%) .  $ p 3$ (Pa) .  Error (%) . 

Exp.  9777.0  X  8700.0  X  299 19.5  X 
C_R  9589.5  −1.9  8982.7  3.2  168 15.4  −43.8 
C_D  9486.1  −3.0  9170.8  5.4  214 76.7  −28.2 
M_D  9523.1  −2.6  9168.5  5.4  247 26.8  −17.4 
F_D  9622.7  −1.6  9212.5  5.9  244 98.2  −18.1 
Case name .  $ p 1$ (Pa) .  Error (%) .  $ p 2$ (Pa) .  Error (%) .  $ p 3$ (Pa) .  Error (%) . 

Exp.  9777.0  X  8700.0  X  299 19.5  X 
C_R  9589.5  −1.9  8982.7  3.2  168 15.4  −43.8 
C_D  9486.1  −3.0  9170.8  5.4  214 76.7  −28.2 
M_D  9523.1  −2.6  9168.5  5.4  247 26.8  −17.4 
F_D  9622.7  −1.6  9212.5  5.9  244 98.2  −18.1 
The time histories of the pressure at the three points for F_D are presented in Fig. 12, as a guide to the discussion on the dynamics of the cavity shedding. The pressures at $ p 1$ and $ p 2$ behave much smoother compared to $ p 3$. All of them show very high pressure when the cavity collapses, corresponding to the minimum bubble volume. This behavior is also observed in the experiments.^{18} Note that this process is violent and may be better approached as a compressible flow, as opposed to the incompressible approach in our current model. As a new cavity begins to form following the collapse of the previous, the pressures stay relatively constant with a slow increase in the total bubble volume. As the time evolves, the pressures begin to drop and the total vapor volume growth rate increases. In addition, the pressure at $ p 3$ fluctuates significantly due to high levels of turbulence in the wake of the cavity. As the vapor volume increases to the maximum, the pressures at $ p 1$ and $ p 2$ drop to the vapor pressure and remain fairly constant until the condensation front passes. The pressures at $ p 1$ and $ p 2$ are very similar as also indicated by the mean values in Table III.
Figure 13 shows the time history of the total pressure at $ p 4$ (location shown in Fig. 14) as the cavity collapses due to massive condensation. Five stages are identified in the process as discussed by Bhatt et al.^{18} in reference to their experiments. The corresponding void fraction contours are shown in Fig. 14. At stage I, $ p 4$ is in the main cavity where the pressure is relatively low. As the condensation front passes $ p 4$ (see Fig. 14), the pressure starts to increase (Fig. 13). There is a sharp decrease in void fraction before and after the condensation front passes $ p 4$, with the pressure also exhibiting a violent transient as indicated in Fig. 13. When the cavity almost disappears at stage III, the pressure is high due to the strong condensation, causing converging flow. As the shedding vortex in the shear layer begins to cavitate again, the void fraction increases in the core of the vortex pushing liquid out with the resulting high pressure at $ p 4$ observed in stage IV. At stage V, the growth of two main cavities (or spanwise vortices) on the two sides of $ p 4$ causes another peak in pressure. In this process, the two spanwise vortices pair with streamwise vortices between them. These results are consistent with the experimental observations, indicating that the model can properly predict the growth, shedding, and condensation processes of the cavity, though the magnitude of the strongest pressure peaks may be affected by the use of an incompressible solver.
D. Shedding frequency
One of the important parameters in this cavitating flow is the cavity shedding frequency. The experiment provides two measurements based on the void fraction and the dynamic pressure. Since the pressure signal is noisier, we choose data based on the void fraction for comparison. The shedding frequency is determined from the frequency with the largest harmonic amplitude in the power spectral density (PSD) of the total vapor volume shown in Fig. 2. Table IV shows the shedding frequencies for all the cases. For the same cavitation number, C_R predicts the lowest value of 15 $ Hz$, while the M_D and F_D cases give the same frequency at 19.5 Hz. (A proper orthogonal decomposition—POD—analysis also presents this as the dominant frequency for the F_D case, see the supplementary material.) As the grid is refined, the frequency seems to converge quickly. To compare with the experimental data, the frequency is expressed in terms of the Strouhal number $ S t = f H / U 0$ and the results are shown in Fig. 15. The overall trend that the frequency decreases for smaller cavitation numbers is predicted by the model. However, the shedding frequency is underestimated, especially for the lower cavitation number. Considering the uncertainty of 0.03–0.04 for the cavitation number reported in experiments, it is possible that the simulated cavitation number is lower.
E. Condensation front
Condensation causes the cavity to collapse, resulting in a moving front with a jump in flow properties, as observed in the experiments and discussed in Bhatt et al..^{18} The $ x \u2212 t$ diagrams of the void fraction at $ y = \u2212 0.25 H$ for F_D are shown in Fig. 16. One thin band of high void fraction, caused by the second cavitation cloud (see also Fig. 28), is formed after the main cloud, as also observed in experiment. Since the width does not increase much in time, the shedding cavity is likely purely transported downstream without significant changes. The main cloud at this vertical position ( $ y = \u2212 0.25 H$) starts to form at around 1H and quickly grows to 6–7H. Two slopes outlined in black correspond to the cavity growth and collapse, indicating that the growth is much faster than the collapse. Figure 17 shows the void fraction and the total pressure at two close times ( $ t = 0.67$ and $ t = 0.68 \u2009 s$) during condensation. There is a pressure jump as the void fraction retreats from $ x = 6 H$ to $ x = 4 H$ as the condensation front advances from downstream to the step. Note that the pressure jump decreases slightly as the condensation front advances. A small peak of 0.1 in the void fraction moves with the condensation front indicating that cavitation occurs after the condensation front; more details are discussed in Sec. V H. The void fraction between $ x = 1 H$ and 4H increases slightly as the condensation front advances, suggesting that there is still certain level of evaporation in the cavity as the pressure is very close to the vapor pressure.
Figure 18 presents the $ x \u2212 t$ diagrams at $ y = \u2212 0.25 H$ for the streamwise velocity, total pressure, and $ d 32$ for one shedding period. As the cavity grows, there are relatively high pressures and streamwise velocities downstream of the cavity in the noncavitating region, due to the evaporation pushing liquid outward. As the condensation front moves toward the step, the velocity there starts to decrease and a pressure discontinuity develops. As the time evolves further, the streamwise velocity becomes negative, pointing toward the step. The pressure jump does not change significantly. Note that the pressure inside the cavity stays close to the vapor pressure as also shown in Fig. 17. The maximum $ d 32$ can exceed 2 mm around 3H and is formed during the cavity collapse. This also indicates that evaporation in some regions inside the cavity takes place even when the overall cavity is collapsing.
In order to obtain the velocity of the condensation front, points around the slope in Fig. 16 with $x$ in the interval $ [ 1 H , \u2009 \u2009 5 H ]$, where the void fraction decreases significantly, are extracted and a linear equation fitting these points is produced. Then, all the slopes for the shedding are averaged to obtain the mean value. As shown in Table V, for the same cavitation number, the differences between methodologies and conditions are small. C_R predicts the highest speed at 0.47 for the same cavitation case. The results are also compared with the experimental data in Fig. 19. The condensation front speed drops as the cavitation number decreases, the same trend shown by the experiments. The overall speed is slightly underestimated for all cases.
Case name .  $ U s / U 0$ . 

C_R  0.47 
C_D  0.42 
M_D  0.40 
F_D  0.41 
$ C _ \sigma L$  0.34 
$ C _ \sigma H$  0.55 
Case name .  $ U s / U 0$ . 

C_R  0.47 
C_D  0.42 
M_D  0.40 
F_D  0.41 
$ C _ \sigma L$  0.34 
$ C _ \sigma H$  0.55 
.  M1 .  M2 .  M3 .  M4 .  M5 .  M6 .  M7 .  M8 .  M9 . 

x/H  2  2  2  3  3  3  4  4  4 
y/H  −0.25  −0.5  −0.75  −0.25  −0.5  −0.75  −0.25  −0.5  −0.75 
.  M1 .  M2 .  M3 .  M4 .  M5 .  M6 .  M7 .  M8 .  M9 . 

x/H  2  2  2  3  3  3  4  4  4 
y/H  −0.25  −0.5  −0.75  −0.25  −0.5  −0.75  −0.25  −0.5  −0.75 
.  M1 .  M2 .  M3 .  M4 .  M5 .  M6 .  M7 .  M8 .  M9 . 

S1  0.470  0.441  0.408  0.433  0.380  0.378  0.441  0.415  0.414 
S2  0.455  0.415  0.401  0.427  0.377  0.370  0.445  0.412  0.397 
.  M1 .  M2 .  M3 .  M4 .  M5 .  M6 .  M7 .  M8 .  M9 . 

S1  0.470  0.441  0.408  0.433  0.380  0.378  0.441  0.415  0.414 
S2  0.455  0.415  0.401  0.427  0.377  0.370  0.445  0.412  0.397 
It is important to note that though the condensation front analyzed at the mixture level (or one fluid level) satisfies the 1D jump condition as a traditional shock does, the mechanism behind them is different. Changes in density of a traditional shock are directly caused by compression resulting from the highpressure forces around the front. For cavitation, the change of mixture density is mainly caused by mild pressure variations around the vapor pressure that result in strong condensation on one side of the front, suggesting that the condensation rate may be critical. 2D incompressible simulations of the condensation front using FLUENT 22R2 with the Zwart–Gerber–Belamri cavitation model^{32} show considerable change in the front speed when the condensation rate is changed. The cavitation model should provide a good estimation of the condensation rate so the front speed compares well to the experimental data for different conditions, as the case with the results shown in Fig. 19.
F. Mach number for the mixture
The mixture density is obtained as $ \rho m = \alpha c \rho c + \alpha d \rho d$ in the following analysis, which essentially equates to $ \alpha c \rho c$ due to the low density of the vapor compared to the liquid. Figure 21 shows the two limiting conditions as predicted by Eq. (23). The speed sound predicted by the homogenous model is very low, and it does not recover the liquid sound speed as the void fraction reaches zero. The ratio of speed of sound with the frozen model with respect to equilibrium model is around 10 to 100 for void fractions between 0.1 and 0.9. This indicates that the speed of sound strongly depends on the model chosen and location in the flow. Figure 22 presents the Mach number based on the local flow velocity (the velocity is equal to the liquid velocity assuming no slip velocity between the two phases) with these two models. It can be seen that the most of the regions result in Mach numbers significantly larger than one when using the homogeneous equilibrium model. On the other hand, the homogeneous frozen model predicts values spanning from the low Mach number to high values larger than one. For the shear region with high void fraction, the Mach number is quite high. In the center of the vortex cavitation, the Mach number is much higher than on the surroundings. If using the velocity at the step, it can be expected that the Mach number will increase significantly near the bottom wall since the velocity there is relatively small.
Budich et al.^{25} estimated the Mach number based on the condensation front velocity for the wedge flow and concluded that it should be larger than one when the front (or shock wave in their paper) forms. We choose the points in Table VI and use the mean flow properties to estimate the Mach number for the two models as shown in Table VIII. The homogeneous equilibrium model predicts a Mach number of around 30, while the homogenous frozen model gives values closer to one. As mentioned by Brennen,^{45} a more realistic situation should lie in between these two models. It is highly possible that the Mach number based on the condensation front velocity is higher than one, consistent with the shock formation condition stated by Budich et al.^{25}
.  M1 .  M2 .  M3 .  M4 .  M5 .  M6 .  M7 .  M8 .  M9 . 

Equilibrium model  39.79  39.17  39.86  32.37  28.76  33.44  37.11  35.56  39.14 
Frozen model  0.96  0.88  0.81  0.89  0.79  0.77  0.89  0.84  0.81 
.  M1 .  M2 .  M3 .  M4 .  M5 .  M6 .  M7 .  M8 .  M9 . 

Equilibrium model  39.79  39.17  39.86  32.37  28.76  33.44  37.11  35.56  39.14 
Frozen model  0.96  0.88  0.81  0.89  0.79  0.77  0.89  0.84  0.81 
G. Vorticity transport equation (VTE)
Figure 23 shows the mean vorticity magnitude in the flow, with a window for averaging indicated by the black rectangle. The norm of each component on the RHS in Eq. (24) is then averaged in the window to investigate the effect to the material derivative of vorticity. Conditions can also be applied to this averaging process. The bottom figure in the panel shows the mean values for the different components on the RHS in Eq. (24). Overall, the vorticity stretching term has the largest contribution. The compressible effect is stronger than the baroclinic effect, as was also observed for the condensation front for the flow over a wedge.^{25} For the buoyancy term, the magnitude is several orders smaller, indicating its negligible effect. The diffusion term also plays an important role due to the strong turbulence produced in the shear layer, which is very different from the inertiadominated flow over the wedge.^{25} For the conditional average over evaporation and high void fraction regions, the contributions for each term are quite similar, indicating that the high void fraction region is related to the evaporation region. In addition, the compressible and baroclinic effects become more important as expected.
The vorticity magnitude during a typical cavity collapse is shown in Fig. 24. The void fraction contour of 0.3 is indicated with black lines. The highest vorticity lies in the turbulent shear region and the reattachment area. As the cavity collapses, the vorticity inside the cavity decreases, which is probably due to the dilation component considering the increasing void fraction in the cavity (see Fig. 29). A clear demarcation occurs at the condensation front. The regions after the front passes have higher vorticity, indicating that the condensation front is a moving source of vortices. The right figures in the panel illustrate the contribution of different terms in VTE for the middle stage on the left. The stretching term has relatively high values behind the step. The compressible and baroclinic terms have large values at the condensation front, supporting that it is a moving source of vorticity. The buoyancy term is quite small in the cavity even though the void fraction is high. The diffusion term is small in the high void fraction regions behind the step but high in the far wake where the turbulence is very strong.
H. Formation of the condensation front
In contrast to reentrant jets, the condensation front forms an advancing jump clearly observable in the mixture flow. Considering the curvature of the flow outside the cavity tail, it is expected that the flow at the reattachment location with high pressure is important for the formation of both processes. The induced bubbly flow may also be sensitive to the flow quantities, as observed by Wu et al.^{22} in their hydrofoil experiments. In the BFS experiments of Bhatt et al., the condensation front is observable as the cavitation number decreases to a certain range. However, they do not provide the flow field when the condensation front forms. Here, we focus on the analysis of the 2D flow as the bubbly front starts to form. A few key time snapshots around the formation of the condensation front are chosen to investigate the mechanisms involved.
Figure 25 presents the flow field in a typical cycle for the front formation. As the cavity grows to its biggest size, almost all the liquid is pushed out of the cavity as indicated by the streamlines, with a recirculation close to the step. No reentrant flow is observed near the bottom wall. When the cavity begins to collapse, small recirculations develop at the cavity tail, while the void fraction inside the cavity increases. Further into the collapse, a clear front forms at the tail of the cavity and the recirculation behind the cavity grows. In the next stage of the cavity collapse, the recirculation keeps growing and eventually breaks down coupled with vortex cavitation in the shear region above. The void fraction inside the remaining cavity continues to increase. This is probably the reason why a reentrant jet cannot form at the bottom wall as is the case of the twist foil,^{10} which is dominated by the reentrant jet. In addition, the void fraction increases slightly near the step, even in the smaller corner recirculation region. This is also consistent with the pressure drop to vapor pressure as indicated by the monitor points discussed in Sec. V C. The front moves toward the step, and the large circulation behind front continues growing. The void fraction inside the remaining cavity increases further, pushing liquid outward. Note that the reattachment location does not change much compared with the initial location.
To quantitatively investigate the near wall flow, a $ y \u2212 t$ diagram of streamwise velocity at a test location of $ x = 3.5 H$ downstream of step is illustrated in Fig. 26. It is clear that the velocity near the wall changes the direction with each cavity shedding cycle, most significantly in the region $ y < 0.4 H$. The time history of the void fraction and streamwise velocity at $ ( x , y ) = ( 3.5 H , \u2009 0.05 H )$ are presented in the right panel of Fig. 26. As the cavity grows and passes the monitor point, the void fraction increases to 0.2–0.3 while exhibiting strong fluctuations. Recalling Eq. (22), the condensation front speed increases as the prefront void fraction decreases, proportional to $ 1 / \alpha d , 1$ for small $ \alpha d , 1$. This indicates that the condensation front speed experiences a large increase if $ \alpha d , 1$ is very small, resulting in a large difference of the condensation front velocity with respect to that in high void fraction regions. This can then destroy the stable front causing instability, preventing the steady moving front. As a result, a high void fraction close to the wall can be a necessary condition for such stable front formation. The streamwise velocity stays fairly constant between 0 and 1 m/s. When the cavity collapses and the condensation front crosses the test location, the void fraction decreases sharply to the background level and the streamwise velocity drops to −3 to −2 m/s. In addition, the velocity fluctuates before the next cavity grows passing the test location due to the highly turbulent flow.
Based on this study, we can conclude that there are important factors influencing the condensation front formation: 1. large pressure at the reattachment location, which initiates the recirculation; 2. high void fraction in the cavity (low cavitation number), causing high void fraction near the bottom wall; 3. considerable evaporation inside the remaining cavity including in the near wall region, which prevents the formation of a reentrant jet flow at the bottom of the cavity. The first two are also discussed for condensation frontdominated flows in the literature.^{11,22,25,27}
I. Cavity evolution history
Figure 27 shows four key times during one shedding cycle for the F_D case. The experimental photographs shown in Fig. 28 along with the simulations are very close in time to the instants shown in Fig. 27, but not exactly. Since the span in experiment is larger than that in simulation, the figures are scaled to provide a better visualization.
As the cavity grows near its maximum rate ( $ t 1$ in Figs. 29 and 30), the low pressure in the vortex shed from the step causes strong evaporation. There are two main evaporation clouds connected with small streamwise cavities. The model predicts well such cavity topology at this particular time, as the vortical structures and isosurfaces of $Q$ clearly show the streamwise vortices between these two clouds. The vapor volume grows to a maximum as the time evolves to $ t 2$, and the main cavity reaches its maximum length and spans over a large space including the region very close to the step. The two main cavity clouds merge into one big cavity. Merging of the two cavities not always happens in the simulations; if the cavities do not merge, the total volume of vapor reaches a lower maximum than if they do, see Fig. 12. At the tail of the cavity, small secondary shedding clouds can be observed in both experiment and simulation, with the presence of vortical structures indicating a highly turbulent flow. Nuclei are attracted to the vortex cores as they react to the vortexinduced pressure gradient, promoting the filaments of cavitation. When the cavity begins to collapse, the main cavity starts shrinking from the tail downstream of the step with a speed slower than the growth, as already discussed. Note that there are a number of small cavities left behind by the condensation front, also present in the experiments. These cavities are caused by strong small vortices (see the right panel in Fig. 28); therefore, a relatively fine grid is necessary to capture them. At time $ t 4$, the vapor volume is at a minimum and the main cavity almost vanished, though there are still a few locations where cavitation occurs despite the fact that the pressure is high in most of the domain; this is also consistent with the experimental observations. Considerable vortical structures are present in the flow, including the spanwise vortices shedding from the step.
Spanwise averaged quantities for the four time steps from Fig. 27 are shown in Fig. 29. The main two cavities formed at $ t 1$ have void fraction around 0.3–0.4. A low void fraction band connecting the cavities is due to cavitation in the core of the streamwise vortices seen in Fig. 28. The total pressures in the two cavities are very low, and the Sauter mean diameter is in the range of 0.6–1.2 mm. At time $ t 2$, the main cavity reaches 60% void fraction in the core with smaller void fraction approaching the step. Some scattered void fraction at the cavity tail can also be observed due to the shedding of cavities produced by random strong vortices. The total pressure is close to vapor pressure in the cavity region increasing sharply at the cavity tail. The Sauter mean diameter continued increasing as the total volume of vapor grew, as a result of bubble coalescence and evaporation in the high void fraction environment. As the cavity shrinks due to condensation (time $ t 3$), a clear condensation front is observed in the void fraction, also present in the experiment. The void fraction in the remaining cavity is very high indicating continued evaporation as mentioned before. Low void fraction locations are left behind by the collapsing main cavity, where strong shear and low total pressure are present. In addition, the high void fraction characterizing the cavity moved closer to the step, consistent with the experimental observations. Though not shown here, a large recirculation (see Fig. 25) formed behind the condensation front and under the shear region. The Sauter mean diameter also grows to 2 mm at the core of the remaining cavity. At time $ t 4$, the void fraction is very low and cannot be discerned with the linear scale in Fig. 29. The total pressure shows high values in most of the domain, but lower pressures are observed in the shear layer. The Sauter mean diameter decreases to the background values in most of the domain.
Figure 30 presents the ratio of DDES length scale $ l \u0303$ to RANS length scale $ l \u0303 RANS$ (see Li et al.^{48} and references therein). A value of 1 means that the region is in the RANS mode and anywhere else is in the LES mode. The highly turbulent areas at and around the shear layer region are in the LES mode, while other parts are in the RANS mode. The near wall regions are in the RANS mode protecting the turbulent boundary layer. Void fraction contour lines shown in Fig. 30 suggest cavitation mostly occurs in LES regions. In particular, the shedding cavity growing at $ t 1$ is in the LES region. Cavitation also occurs in near wall RANS regions. In addition, the region close to the step is in the RANS mode since no strong turbulence is produced there, where cavitation occurs when the cavity collapses at time $ t 3$.
VI. SUMMARY AND CONCLUSIONS
The incompressible polydisperse twophase flow model of Li and Carrica^{10} is used to study the cavitating flow over a backward facing step with RANS and DDES. A few important model parameters are investigated to study their influence on predictions. Overall, the cavitation model can predict the partial cavitation regime investigated, with satisfactory results for void fraction, shedding frequency, condensation front speed, and the total pressure in the flow field. The RANS model can predict the unsteady shedding behavior, but overpredicts the cavity length and void fraction. For the maximum void fraction location, shedding frequency, and condensation speed, differences are not significant between the models for the same cavitation number. The total pressure at locations close to the step is also model independent, and predictions agree well with the experimental data. The total pressure close to the reattachment location is more difficult to predict and varies significantly between simulated cases. In general, DDES presents better results than RANS and finer grids improve the predictions, resolving the smaller cavitation structures induced by strong vortices. The model also predicts the same trends as the experiments in terms of the cavity length, shedding frequency, and condensation front as the cavitation number changes, though the results show considerable differences with respect to the experimental data for some quantities. Overall, the model can capture well the main flow features without modeling constants for the cavitation source.
The condensation front velocity, approximated by a 1D analytic equation as done in Bhatt et al.,^{18} is consistent with the experiments. The cavitation model, though incompressible for each fluid, predicts the occurrence of condensation fronts for the different cavitation numbers explored and under different numerical and modeling conditions, indicating that the formation of this front is not associated with the traditional compressible nature of the mixture and the pressure wave generated by the cavity collapse, but with the changes of mixture density caused by evaporation and condensation. The formation of the condensation front is driven by strong condensation rates, suggesting the importance of a proper condensation rate model when studying transient cavitating flows. In addition, the results suggest that the use of a compressible solver for the flow is not a necessary condition to simulate the formation of a condensation front.
Two extreme models to compute the sound speed in the mixture as a postprocessing step are investigated. The resulting Mach number using the local flow velocity shows that the homogeneous equilibrium model predicts significantly higher values than those by the homogeneous frozen model, consistent with results shown with the experiments. The estimated local Mach number, based on the condensation front velocity, is shown to be close to one for the homogenous frozen model.
The 2D flow analysis for cavity collapse shows that considerable evaporation and high void fraction exist in the remaining of the cavity, including near the bottom wall region. The POD analysis for the 2D flow, presented as the supplementary material, shows that the first few modes dominate the flow and multiple frequencies are observed in the modes, but it is dominated by the cavity shedding frequency. The vorticity transport equation for the liquid is analyzed numerically including the buoyancy term. It is found that the condensation front is a moving source of the vorticity. The “compressible” (in reality due to changes in void fraction) and baroclinic effects are considerable near the jump, while the buoyancy term is several orders of magnitude smaller and negligible. The vorticity inside the remaining cavity is found to decrease as the cavity collapses.
The polydisperse cavitation model has the capability to compute the bubble size distribution, a unique advantage to this approach. The results show large differences in bubble sizes for different locations, indicating that the evaporation/condensation source can significantly differ in regions with similar pressures because the source is strongly dependent on the bubble size. Smaller bubbles, resulting from bubble fission during the cavity collapse, are located mostly in the downstream region of the cavity, where collapse is more violent. Larger bubbles reside mostly closer to the step, with the largest residing nearest the step, where the flow is milder. While qualitative comparison shows consistency with videos of the experiments for the locations of larger and smaller bubbles, the lack of experimental bubble size distributions for the validation results in large uncertainty for the size predictions.
SUPPLEMENTARY MATERIAL
See the supplementary material for details of an animation of the DDES simulation with the fine grid (video I) and a section discussing results of a POD analysis.
ACKNOWLEDGMENTS
This research was sponsored by the U.S. Office of Naval Research through MURI Grant No. N00014172676, University of Minnesota lead institution, Dr. KiHan Kim and Julie Young program officers.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Jiajia Li: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Pablo M. Carrica: Conceptualization (equal); Funding acquisition (lead); Methodology (equal); Project administration (lead); Resources (lead); Software (equal); Supervision (lead); Writing – review & editing (lead).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.