Comparative performance assessment between incompressible and compressible solvers to simulate a cavitating wake

Citation: Chen, J.; Geng, L.; Jou, E.; Escaler, X. Comparative Performance Assessment between Incompressible and Compressible Solvers to Simulate a Cavitating Wake. Fluids 2024,9, 218. https://doi.org/10.3390/fluids9090218 Academic Editors: Andrei Lipatnikov and Manolis Gavaises Received: 31 July 2024 Revised: 28 August 2024 Accepted: 10 September 2024 Published: 18 September 2024 Copyright: © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). fluids Article Comparative Performance Assessment between Incompressible and Compressible Solvers to Simulate a Cavitating Wake Jian Chen 1, Linlin Geng 2, Esteve Jou 1and Xavier Escaler 1,* 1Barcelona Fluids & Energy Lab, Universitat Politècnica de Catalunya, 08028 Barcelona, Spain; [email protected] (J.C.); [email protected] (E.J.) 2National Research Center of Pumps, Jiangsu University, Zhenjiang 212013, China; [email protected] *Correspondence: xavier[email protected] Abstract: To study the effects of fluid compressibility on the dynamics of a cavitating vortex street flow in a regime where the vortex shedding frequency increases as a result of the cavitation increase, the cavitating wake behind a wedge was simulated employing both incompressible and compressible solvers. To do this, a compressible cavitation model was implemented, modifying the Zwart-Gerber-Belamri (ZGB) incompressible solver and including a pressure limit and absorbing boundary conditions to prevent a non-physical pressure field. To validate the performance of the compressible model, preliminary simulations were carried out on a 1D Sod cavitating tube and the cavitating vortex shedding behind a circular body at laminar flow conditions. The results of the cavitating wake behind the wedge with the incompressible and the compressible solvers showed similar predictions in terms of pressure, vortex shedding frequency, and instantaneous and average vapor volume fraction profiles. In spite of this, differences were obtained in the energy content of the fluid force fluctuations on the body at higher frequencies, which appear to be better resolved and amplified when the compressibility model is considered. Overall, both solvers provided comparable results in terms of cavitation phenomena that are well aligned with experimental observations. Keywords: vortex street; cavitation; fluid compressibility; sponge layer; spectral content 1. Introduction Vortical flow structures and cavitation are two common phenomena encountered in hydraulic machinery [ 1 – 3 ]. As the fluid passes the fixed vanes, the guide vanes, or the runner blades, it creates complex wake flow patterns that might develop in what is known as a vortex street. The resulting alternating shed vortices induce force fluctuations on the solids and provoke vortex-induced vibrations that may induce severe fatigue damage. In addition, vaporous cavities will be formed and persist inside the center of the alternating shed vortices if the ambient pressure is low enough [ 4 ]. As these cavities move to the high-pressure region and collapse, they can also induce severe material erosion on adjacent surfaces. It is well known that the occurrence and development of cavitation can strongly alter the dynamic behavior of vortex street flows [ 5 – 8 ]. For example, it was reported that cavitation tends to enhance the amplitude of vortex-induced vibrations [ 7 ]. Therefore, further exploration of the effect of cavitation on the vortex street flow is essential to mitigate the damage caused by the cavitation and/or vortex-induced vibrations. Currently, the computational fluid dynamic (CFD) method has been proved to provide detailed insights into complex flow fields, particularly in cases involving cavitation. Various numerical approaches, each with different levels of complexity, are available to simulate cavitating flow. For a comprehensive review of cavitation modeling approaches, see Nied´zwiedzka et al. [ 9 ] and Folden and Aschmoneit [ 10 ]. Among these numerical models, the transport equation model (TEM) is one of the most widely used; this includes an additional transport equation for the vapor phase volume fraction and a source term that Fluids 2024,9, 218. https://doi.org/10.3390/fluids9090218 https://www.mdpi.com/journal/fluids Fluids 2024,9, 218 2 of 20 accounts for the mass transfer between phases. Without considering fluid compressibility, TEM-based simulations can effectively capture the transition from sheet to cloud cavitation, driven by a liquid re-entrant jet across different geometries [ 11 – 13 ], aligning well with experimental observations [ 1 , 14 – 17 ]. Conversely, the transition from sheet cavitation to cloud cavitation can be dominated by condensation shock [ 18 – 20 ], which is highly related to the fluid compressibility of the mixture. Coupled with the equation of state for the liquid and vapor, TEM can successfully capture the condensation shock and its impact on the transition from sheet cavitation to cloud cavitation [13,20–23]. For the developed cavitation within the vortex street flows, bubble–vortex interactions alter the growth and collapse dynamics of the cavitation [ 24 – 30 ]. This makes the roles of fluid compressibility in cavitating vortex street flows distinct from those observed in unsteady cloud cavitation. Previous numerical investigations highlighted the impact of fluid compressibility on cavitating vortex street flows, particularly at extremely low cavitation numbers. For example, Brandao and Mahesh [ 31 ] numerically captured the condensation shock on the surface of the circular cylinder where cavitation was highly developed. Additionally, Kim et al. [ 32 ] confirmed that fluid compressibility can alter the morphology of fully developed cavities when the cavitation number is less than half of the cavitation inception number. Recently, Wang et al. numerically captured the interaction between the microbubbles and streamwise vortices in cavitating vortex street flow using the Euler–Lagrange method [ 33 ]. However, the effects of fluid compressibility on the dynamics of vortex street flows have received limited attention, especially at intermediate cavitation numbers. The aim of this study is to demonstrate the influence of fluid compressibility on the dynamics of a cavitating vortex street by numerically resolving the cavitating flow field using both incompressible and compressible cavitation solvers. This paper is composed of five sections. Following this introduction, the second section explains the governing equations. The third section details the corresponding numerical treatment. In the fourth section, the implemented numerical solvers are validated, and the results from both the incompressible and compressible solvers are compared. Finally, the fifth section presents the conclusions and summarizes the results. 2. Governing Equations 2.1. Mass Conservation Equation The equation for the conserved mass without an external mass source is given by: ∂ρ/∂t+∇·(ρu)=0 (1) where ρ is the density, u is the velocity, and t is the time. This equation is valid for both incompressible and compressible flows. 2.2. Momentum Conservation Equation In an inertial reference frame, the conserved momentum equation can be written as follows: ∂(ρu)/∂t+∇·ρuOu=−∇p+∇·τ(2) where pis the pressure, and τis the stress tensor, which is defined as: τ=µh∇u+∇uT−2/3(∇·u)Ii(3) where µis the fluid molecular viscosity, and Iis the unit tensor. Fluids 2024,9, 218 3 of 20 2.3. Realizable k-Epsilon Delayed Detached Eddy Simulation (DDES) Model In the realizable k -epsilon DDES model, the governing equation for the turbulence kinetic energy, k, can be written as follows: ∂(ρk)/∂t+∇·(ρuk)=Pk−Y* k+∇·[(µ+µt/σk)∇k](4) where ε is the dissipation, Pk is the generation term of k , µt is the turbulent viscosity, and σk is the turbulent Prandtl number. Furthermore, the dissipation of turbulent kinetic energy, Y* k, is given by: Y* k=ρk3 2 ldes (5) where ldes =min(lrke,lles);lrke =k3 2 ε;lles =Cdes∆max (6) where lrke and lles are the calculated length scales based on the realizable k -epsilon model and LES, respectively. ∆max is the maximum length of the cell edge, and Cdes = 0.61 is the model constant. If lrke =lles,ldes is replaced by the following functions: ldes =lrke −fymax(0, lrke −Cdes∆max)(7) fy=1−tan hh20.0ry3i(8) ry=(µ+µt)/ρ √u:uκy2(9) where yrepresents the wall distance, and κ= 0.41. 2.4. Cavitation Modeling The volume fraction transport equation for the vapor volume fraction, αv , is given by: ∂(ρvαv)/∂t+∇·(ρvαvu)=. m(10) The ρand µof the fluid mixture are defined, respectively, as: ρ=αvρv+(1−αv)ρl(11) µ=αvµv+(1−αv)µl(12) where ρv and ρl are the water and vapor densities, respectively, and µv and µl are the water and vapor dynamic viscosities, respectively. The mass transfer rate between the liquid and vapor due to the cavitation is modeled with a source term, . m , presented in Equation (10). If the effects of viscosity, non-condensable gas, surface tension, and second-order derivation are neglected, the Rayleigh–Plesset equation can be simplified and written as: . R=v u u t2Pre f −pv 3ρl (13) where pis the saturated vapor pressure, and Ris the bubble radius. Assuming that the distribution of the nucleation inside the flow field is uniform and independent of the flow structures, then the relationship between the bubble radius, R , and αvis given by: αv=nv4π 3R3(14) Fluids 2024,9, 218 4 of 20 where nvis the number of nucleations per unit volume. Then, . mis calculated by: . m=ρv . αv=nvρv34π 3R2. R(15) By submitting n given by Equations (13) and (14) into Equation (15), . m can be expressed as: . m=3ρvαv Rv u u t2Pre f −pv 3ρl (16) In the vaporization process, the nucleation site density must decrease as the vapor volume fraction increases. To model this process, α needs to be replaced with (1−αv)αnuc in Equation (16). Then, the final ZGB cavitation model [34] is presented as follows: . m=     Fc3αvρv R0q2 3 (p−pv) ρlp>pv Fv3ρv(1−αv)αnuc R0q2 3 (pv−p) ρlp<pv (17) where pv is the saturated vapor pressure. The initial value of the bubble radius is R0=1µm, and the nucleation site of the volume fraction αnuc = 5 × 10 −4 . Here, the optimal empirical condensation and vaporization coefficients have been selected as Fc= 0.001 and Fv=50.0, respectively. 2.5. Equation of State The density of the compressible liquid phase can be computed using the Tait equation [22,35]: ρ(p)l=ρl,satp+B pv+B1/N (18) and the corresponding sound speed is defined as c(p)l=s∂p ∂ρ =sN ρl,sat (pv+B)p+B pv+B(N−1)/N (19) where B is the water bulk modulus, 3.1 × 10 8 Pa, ρl,sat is the liquid density at pv , 998.18 kg/m3, and Nis a model constant equal to 7.15. The density of the compressible vapor phase using the polytropic equation of state [ 36 ] can be written as: ρ(p)v=(p/Cv)1/n(20) and the corresponding sound speed is: c(p)v=s∂p ∂ρ =qCvn(p/Cv)(n−1)/n(21) where the model constant Cv can be estimated from a density of 0.017 kg/m 3 at pv . The polytropic exponent nfor an adiabatic assumption is 1.4. 2.6. Sponge Layer Conditions To date, although the use of the pressure-based method has achieved notable success in a wide range of compressible cavitating flows [ 37 , 38 ], several unresolved numerical difficulties remain, particularly concerning the determination of the inlet and outlet boundary conditions in CFD simulations. Typically, such calculations are conducted on a truncated domain of the entire system, and the standard inlet/outlet boundary conditions may fail to allow some of the flow features to leave the computational domain as physically expected. Fluids 2024,9, 218 5 of 20 For that reason, artificial treatments of the truncated domains are required to avoid pressure waves reflected at the boundaries. Therefore, the sponge layer was adopted for its flexibility and simplicity [ 39 ]. The set of formulas of the sponge layer, based on the pressure and velocity for the pressure-based method, is given as follows: ∂ρ/∂t+∇·(ρu)=σ(x)∂ρ ∂pp=pre f pre f −p(22) ∂(ρu)/∂t+∇·(ρuu +pI−τ)=σ(x)ρure f −u(23) where the subtitle “ref” refers to a target value of the flow variable, and σ(x) is the function of the damping coefficient, which is defined as: σ(x)=3αpol(x/L)βpol 2∆tγpol (24) where αpol = 1.0, βpol = 3.0, γpol = 1.0, xis the distance from the leading edge of the sponge layer, and Lis the length of the sponge layer, as shown in Figure 1. Fluids 2024, 9, x FOR PEER REVIEW 5 of 21 2.6. Sponge Layer Conditions To date, although the use of the pressure-based method has achieved notable success in a wide range of compressible cavitating flows [37,38], several unresolved numerical difficulties remain, particularly concerning the determination of the inlet and outlet boundary conditions in CFD simulations. Typically, such calculations are conducted on a truncated domain of the entire system, and the standard inlet/outlet boundary conditions may fail to allow some of the flow features to leave the computational domain as physically expected. For that reason, artificial treatments of the truncated domains are required to avoid pressure waves reflected at the boundaries. Therefore, the sponge layer was adopted for its flexibility and simplicity [39]. The set of formulas of the sponge layer, based on the pressure and velocity for the pressure-based method, is given as follows: 𝜕𝜌𝜕𝑡 ⁄+∇∙󰇛𝜌𝒖󰇜=𝝈󰇛𝒙󰇜𝜕𝜌 𝜕𝑝𝑝−𝑝 (22) 𝜕󰇛𝜌𝒖󰇜𝜕𝑡 ⁄+∇∙󰇛𝜌𝒖𝒖+𝑝𝑰−𝝉󰇜=𝝈󰇛𝒙󰇜𝜌𝒖−𝒖 (23) where the subtitle “ref” refers to a target value of the flow variable, and 𝝈󰇛𝒙󰇜 is the function of the damping coefficient, which is defined as: 𝝈󰇛𝒙󰇜=3𝛼󰇛𝒙𝐿 ⁄󰇜 2∆𝑡 (24) where 𝛼=1.0 , 𝛽=3.0 , 𝛾=1.0, 𝑥 is the distance from the leading edge of the sponge layer, and 𝐿 is the length of the sponge layer, as shown in Figure 1. Figure 1. Schematic of the sponge layer implementation. 3. Numerical Method In this section, the governing equations are discretized using the finite volume method (FVM) over the collocated unstructured mesh. More attention will be paid to the numerical algorithm for the pressure-based multiphase solver with/without consideration of fluid compressibility using a predefined macro inside the ANSYS Fluent 2022R2 platform provided as a Supplementary File. 3.1. Incompressible Mixture/VOF Model Without consideration of fluid compressibility, the fluid densities are constant. The numerical algorithm for the mass conservation and volume fraction transport equation with linearized source terms is detailed in the following subsection. Figure 1. Schematic of the sponge layer implementation. 3. Numerical Method In this section, the governing equations are discretized using the finite volume method (FVM) over the collocated unstructured mesh. More attention will be paid to the numerical algorithm for the pressure-based multiphase solver with/without consideration of fluid compressibility using a predefined macro inside the ANSYS Fluent 2022R2 platform provided as a Supplementary File. 3.1. Incompressible Mixture/VOF Model Without consideration of fluid compressibility, the fluid densities are constant. The numerical algorithm for the mass conservation and volume fraction transport equation with linearized source terms is detailed in the following subsection. 3.1.1. Incompressible Volume Continuity Equation To guarantee mass conservation and numerical stability, the pressure-correction equation is based on the total volume continuity instead of the mass conservation equation. 1 ρv∂(ρvαv)/∂t+∇·(ρvαvu)−. m+1 ρl∂(ρlαl)/∂t+∇·(ρlαlu)+. m=0 (25) Integrating Equation (25) over a control volume, the discretized form of the volume continuity is given by: Fluids 2024,9, 218 6 of 20            1 ρv ρvαv−ρvαv0 ∆tdV +∑ f ρvαvufAf!+ 1 ρv ρlαl−ρlαl0 ∆tdV +∑ f ρlαlufAf!           =                1 ρv−1 ρl". m*−∂. m ∂p*p*−pv#+ ∂. m ∂p*1 ρv−1 ρl | {z } ds/dp (p−pv)                (26) where dV is the volume of the control cell, Af is the face area vector within the control volume, and ρvαvuf and ρlαluf are the vapor and liquid mass fluxes, respectively. On the collocated grid, the phase mass flux is computed using a Rhie–Chow interpolation. More importantly, the term ds/dp is the linearized mass transfer related to the pressure for the pressure correction equation, which can enhance the numerical stability of the cavitation simulation. 3.1.2. Incompressible Second Phase Fraction Equation In ANSYS Fluent, only the transport equation for the second phase volume fraction is solved using the VOF or mixture model. Normally, the vapor phase is set to be the second phase. Thus, the vapor volume fraction is obtained with the vapor phase continuity equation, as presented in Equation (10). Correspondingly, the discretized form of the equation can be written as follows: 1 ρv ρvαv−ρvαv0 ∆tdV +∑ f ρvαvufAf!= . m ρv dV (27) Assuming that the mass source term . mcan be rewritten as the function of αvand αl: . m ρv =. m0+Slαl−Svαv(28) With the linearized mass source term . m, the Equation (27) can be rewritten as: 1 ρv ρvαv−ρvαv0 ∆tdV +∑ f ρvαvufAf!=Sc+SpαvdV (29) where Sc=. m0+Sl Sp=−(Sv+Sl)(30) where Sp is the linear part of the source term . m , and Sc is the part of . m that cannot be linearized. 3.2. Compressible Mixture/VOF Model In this section, the robust numerical algorithm for multiphase flows to handle fluid compressibility is introduced [ 37 ]. With the SIMPLE method, velocity and density fields can be decomposed into two parts, as follows:        ρv=ρv*+ρv′=ρv*+∂ρv* ∂pp′ ρl=ρl *+ρl′=ρl *+∂ρl * ∂pp′ u=u*+u′ (31) where “*” indicates the tentative values, and “′” indicates the variable correction. Fluids 2024,9, 218 7 of 20 3.2.1. Compressible Volume Continuity Equation Unlike the incompressible pressure-correction equation, the phase masses for vapor and liquid, ρvαv and ρlαl , respectively, and the phase flux terms for vapor and liquid, ρvαvu and ρlαlu, respectively, are computed with the following equations:    ρvαv=ρv*αv+∂ρv* ∂pp′αv ρlαl=ρl*αl+∂ρl* ∂pp′αl (32)                        ρvαvu=ρv*+ρv′αvu*+u′=       ρv*u*αv+∂ρv* ∂pp′u*αv+ ρv*u′αv+ρv′αvu′ | {z } high−order        ρlαlu=ρl*+ρl′αlu*+u′=       ρl*u*αl+∂ρl* ∂pp′u*αl+ ρl*u′αl+ρl′αlu′ |{z} high−order        (33) Normally, the high-order correction terms in Equation (33) are omitted due to their small magnitude compared with the rest of the terms and their negligible influence on the resolved field. Then, the discretized form of the volume continuity can be written as:                                    1 ρv              ρv*αv*+∂ρv* ∂pp′−ρv0αv0 ∆tdV+ ∑ f ρv,f*u*αv+∂ρv* ∂pf p′u*αv+ρv,f*u′αv!Af              + 1 ρl              ρl*αl*+∂ρl* ∂pp′−ρl0αl0 ∆tdV+ ∑ f ρl,f*u*αl+∂ρl* ∂pf p′u*αl+ρl,f*u′αl!Af                                                 =                1 ρv−1 ρl". m*−∂. m ∂p*p*−pv#+ +∂. m ∂p*1 ρv−1 ρl | {z } ds/dp (p−pv)                (34) where ρv,f *u*αvis the phase mass flux at the cell face. The above equation not only couples the pressure and velocity but also involves the interaction between the flow field and the mass transfer source term. Note that the terms ∂ρv* ∂p and ∂ρl* ∂p are related to the effect of the vapor and liquid compressibility on the pressure correction, respectively. By default, these terms are approximated by the upwind scheme. The compressibility-related term in the transient term and mass transfer source term can be easily linearized, which can significantly enhance the numerical stability without affecting the final solution. 3.2.2. Compressible Second Phase Fraction Equation With consideration of the compressible effect of the second phase, the discretized form of the volume fraction equation can be written as follows when the second phase is the vapor: 1 ρv ρvαv−ρv0αv0 ∆tdV +∑ f ρv,fαvufAf!=Sc+SpαvdV (35) And the left-hand side of the equation can be rewritten as: Fluids 2024,9, 218 8 of 20 αv−αv0 ∆tdV +∑ f αvufAf=                  Sc+SpαvdV+              −ρv−ρv0 ρv αv0dV ∆t+ ∑ f ρv−ρv,f ρv αvufAf! |{z } expansion source                               (36) where the last term is the expansion source term accounting for the effect of fluid compressibility on the volume conservation equation. It is composed of two parts: one corresponding to the transient term and the other one to the advection term. 3.3. Pressure Limits As we know, the pressure-based solver can provoke an unbounded pressure field, i.e., a negative pressure field. For the compressible flow, it is necessary to limit the pressure since the obtained negative pressure field is out of the predefined range for the fluid material. However, limiting the pressure is a challenging numerical issue, which will lead to a decrease in the convergence rate and to numerical instability. In the case of multiphase compressible flows, the situation will be more severe. As mentioned by Li and Vasquez [ 37 ], it is difficult to directly limit the pressure above zero in the region where only a negligible amount of gas is present. In such a case, the limitation of the pressure must consider the effect of the local flow characteristic. According to the pressure-limited method proposed by Li and Vasquez [ 37 ], the fluid density is computed using the barotropic law, with a pre-described pressure limit when the local pressure turns negative. For instance, the Tait equation for the compressible liquid density with the pressure limited can be rewritten as: ρ(p)l=ρl,satmax(p,plim)+B pv+B1/N (37) and the corresponding sound speed with the pressure limited is defined as c(p)l=s∂p ∂ρ =sN ρl,sat (pv+B)max(p,plim)+B pv+B(N−1)/N (38) where plim is the pre-described pressure limit. 4. Validation 4.1. Case 1: 1-D Two-Phase Time-Dependent Test Case In this section, the compressible cavitation model is used to describe the effects of a density reduction below the density value of the saturated liquid induced by two symmetrical expansion waves, as in the numerical setup defined by Schmidt [ 40 ]. The domain consists of a 1-D tube with a length of 1 m. At time t = 0 s, the tube is full of pure liquid water, and the pressure is constant with a value of p = 0.9 bar. The velocity field is assumed to jump at x = 0.5 m from the left velocity, uL = − 10 m/s, to the right velocity, uR= 10 m/s, to force the phase change from liquid to vapor to occur. Furthermore, the domain is divided into 1000 cells (four equally spaced cells in the y direction and two-hundred-fifty in the x direction), and the time integration is performed using the first-order implicit scheme with a time step ∆t= 1 ×10−7s. To verify the implementation of the compressible cavitation model, the current numerical results were compared with the ones obtained by Schmidt [ 40 ]. Figure 2presents the pressure, the velocity, the vapor volume fraction, and the sound speed distributions along the tube at time t = 1.5 × 10 −4 s using the present compressible cavitation model and the results obtained by Schmidt [ 40 ]. Without experiencing numerical oscillations in regions with high gradients, the currently implemented solver is proved to be numerically stable. Furthermore, the distribution of the flow quantities along the tube is almost identical Fluids 2024,9, 218 9 of 20 between the current ones and those extracted from the reference simulation. The precise agreement between both sets of results demonstrates the validity of the implemented compressible cavitation model. Fluids 2024, 9, x FOR PEER REVIEW 9 of 21 To verify the implementation of the compressible cavitation model, the current numerical results were compared with the ones obtained by Schmidt [40]. Figure 2 presents the pressure, the velocity, the vapor volume fraction, and the sound speed distributions along the tube at time 𝑡 = 1.5 × 10−4 s using the present compressible cavitation model and the results obtained by Schmidt [40]. Without experiencing numerical oscillations in regions with high gradients, the currently implemented solver is proved to be numerically stable. Furthermore, the distribution of the flow quantities along the tube is almost identical between the current ones and those extracted from the reference simulation. The precise agreement between both sets of results demonstrates the validity of the implemented compressible cavitation model. Figure 2. Comparison of flow quantities obtained with the present model and the simulation by Schmidt (2015) [40] at time 𝑡 = 1.5 × 10−4 s for the 1-D two-phase time-dependent case. 4.2. Case 2: Cavitating Flow over a Circular Cylinder The present simulations of the non-cavitating and cavitating flows over a stationary cylinder were performed with a large 2D computational domain of dimensionless dimensions relative to the cylinder diameter, D, in the horizontal and vertical ranges −50 ≤ 𝑥/𝐷 ≤ 50 and −50 ≤ 𝑦/𝐷 ≤ 50, respectively, to avoid the effect of the domain boundaries, as depicted in Figure 3. The cylinder is located in position (0, 0), and the number of grid elements in the circumferential and radial directions are 480 and 220, respectively. The symbol 𝜃 represents the angle around the cylinder surface measured from its front stagnation point. The inflow condition of 𝑢=𝑢, and the pressure condition of 𝑝 = 𝑝 is set at the boundaries located at 𝑥/𝐷 = −50 and 𝑥/𝐷 = 50, respectively, while the symmetry condition is used at the top and bottom boundaries. Figure 2. Comparison of flow quantities obtained with the present model and the simulation by Schmidt (2015) [40] at time t= 1.5 ×10−4s for the 1-D two-phase time-dependent case. 4.2. Case 2: Cavitating Flow over a Circular Cylinder The present simulations of the non-cavitating and cavitating flows over a stationary cylinder were performed with a large 2D computational domain of dimensionless dimensions relative to the cylinder diameter, D, in the horizontal and vertical ranges −50 ≤x/D≤50 and − 50 ≤y/D≤ 50, respectively, to avoid the effect of the domain boundaries, as depicted in Figure 3. The cylinder is located in position (0, 0), and the number of grid elements in the circumferential and radial directions are 480 and 220, respectively. The symbol θ represents the angle around the cylinder surface measured from its front stagnation point. The inflow condition of u=ure f , and the pressure condition of p=pre f is set at the boundaries located at x/D = − 50 and x/D = 50, respectively, while the symmetry condition is used at the top and bottom boundaries. The drag and lift coefficients, denoted as CDand CL, are defined as: CD=Fx 1 2ρlUre f 2D(39) CL=Fy 1 2ρlre f 2D(40) where Fx and Fy are the streamwise and transverse components of the force acting on the cylinder surface, respectively, and Ure f is the free-stream velocity. The unsteady vortex shedding case is calculated for liquid water at Re = 200 as in previous studies [ 21 , 41 – 47 ]. Vortices are generated on the cylinder surface and detach alternatively from its upper and lower surfaces, which results in fluctuations of CD and CL as shown in Figure 4. The time average value of CD , denoted as CD,av , the maximum value of CL , denoted as CL,max , and the St value are listed in Table 1and compared with Fluids 2024,9, 218 16 of 20 Fluids 2024, 9, x FOR PEER REVIEW 16 of 21 (a) (b) Figure 13. (a) Lift coefficient time histories and (b) corresponding spectra with the incompressible and the compressible cavitation models at 𝜎 = 1.9. 5.3.3. Cavitation Structures Figure 14 shows a visualization of the 3D cavitation structures inside the shed vortices using iso-surface plots, with the numerical results obtained with both the incompressible and the compressible cavitation solvers. They are compared in the same figure with the photographs taken during the experimental tests carried out by Wu et al. [8]. It can be seen that the first few pairs of spanwise vortex cores (identified with the 𝑄 criterion) are filled with vapor, and these cavity structures (represented with iso-surfaces of 𝛼 = 0.05) are advected downstream. Furthermore, the relative positions of the vortices and the spacings obtained with the incompressible and compressible cavitation models are quite similar, and they align well with the experimental results. (a) (b) (c) Figure 14. Comparisons of (a) the experimentally [8] and the numerically obtained cavity structures using (b) the incompressible and (c) the compressible cavitation models at 𝜎 = 1.9. The vortical structures are identified with an iso-surface of 𝛼 = 0.05, and the vorticity level is indicated with the 𝑄 criterion. Figure 14. Comparisons of (a) the experimentally [ 8 ] and the numerically obtained cavity structures using (b) the incompressible and (c) the compressible cavitation models at σ = 1.9. The vortical structures are identified with an iso-surface of αv = 0.05, and the vorticity level is indicated with the Qcriterion. Figure 15 shows the simulated instantaneous field of αv in the cavitating wake at different instants of time obtained with the incompressible (Incomp.) and compressible (Comp.) cavitation solvers at σ = 1.9 as well as the experimental (Exp.) images reported by Wu et al. [ 8 ] using time-resolved X-ray densitometry. All of them show clearly the periodically shed vortices at the vortex shedding frequency. The evolution of the cavitation structures with time confirms that the periodic shedding is governed by the alternating shedding vortex street behind the wedge. In each shedding period, a cavity starts to form and fill the center of the shed vortex. Then, the cavitating vortex is advected from the attached boundary layer and shed from the trailing edge. The alternating cavitating shedding vortices form the vortex street behind the wedge. The comparison between numerical results using incompressible and compressible cavitation solvers and the experimental results indicates that both cavitation solvers can capture the main flow characteristics in terms of the cavities morphology and their relative position, as can be seen in Figure 15a,c. Figure 16 shows the time average void fraction field calculated with the incompressible (Incomp.) and compressible (Comp.) cavitation solvers and the experimental results. It can be seen that the simulations with two different cavitation solvers produce almost indistinguishable mean and RMS void fraction fields behind the wedge. The shape of mean cavity structures obtained numerically is comparable to the pattern observed in the experimental image. Note that there is a downstream offset in the mean and RMS void fraction fields between the numerical and the experimental results, which may be due to the low numerical resolution of the locally unstable flow structures within the detachment area. Fluids 2024,9, 218 17 of 20 Fluids 2024, 9, x FOR PEER REVIEW 17 of 21 Figure 15 shows the simulated instantaneous field of 𝛼 in the cavitating wake at different instants of time obtained with the incompressible (Incomp.) and compressible (Comp.) cavitation solvers at 𝜎 = 1.9 as well as the experimental (Exp.) images reported by Wu et al. [8] using time-resolved X-ray densitometry. All of them show clearly the periodically shed vortices at the vortex shedding frequency. The evolution of the cavitation structures with time confirms that the periodic shedding is governed by the alternating shedding vortex street behind the wedge. In each shedding period, a cavity starts to form and fill the center of the shed vortex. Then, the cavitating vortex is advected from the attached boundary layer and shed from the trailing edge. The alternating cavitating shedding vortices form the vortex street behind the wedge. The comparison between numerical results using incompressible and compressible cavitation solvers and the experimental results indicates that both cavitation solvers can capture the main flow characteristics in terms of the cavities morphology and their relative position, as can be seen in Figure 15a,c. (a) 𝑡 (b) 𝑡 = 𝑡 + 3 ms (c) 𝑡 = 𝑡 + 5 ms (d) 𝑡 = 𝑡 + 8 ms Fluids 2024, 9, x FOR PEER REVIEW 18 of 21 (e) 𝑡 = 𝑡 + 10 ms Figure 15. Comparisons of the experimentally (Exp.) [8] and numerically obtained void fraction contour plots at different instants of time during the vortex shedding using the incompressible (Incomp.) and the compressible (Comp.) cavitation models at 𝜎 = 1.9. Figure 16 shows the time average void fraction field calculated with the incompressible (Incomp.) and compressible (Comp.) cavitation solvers and the experimental results. It can be seen that the simulations with two different cavitation solvers produce almost indistinguishable mean and RMS void fraction fields behind the wedge. The shape of mean cavity structures obtained numerically is comparable to the pattern observed in the experimental image. Note that there is a downstream offset in the mean and RMS void fraction fields between the numerical and the experimental results, which may be due to the low numerical resolution of the locally unstable flow structures within the detachment area. (a) (b) Figure 16. Comparisons of the experimentally (Exp.) [8] and numerically obtained time average (a) and (b) RMS void fraction fields using the incompressible (Incomp.) and the compressible (Comp.) cavitation models at 𝜎 = 1.9. Figure 15. Comparisons of the experimentally (Exp.) [ 8 ] and numerically obtained void fraction contour plots at different instants of time during the vortex shedding using the incompressible (Incomp.) and the compressible (Comp.) cavitation models at σ= 1.9. Fluids 2024,9, 218 18 of 20 Fluids 2024, 9, x FOR PEER REVIEW 18 of 21 (e) 𝑡 = 𝑡 + 10 ms Figure 15. Comparisons of the experimentally (Exp.) [8] and numerically obtained void fraction contour plots at different instants of time during the vortex shedding using the incompressible (Incomp.) and the compressible (Comp.) cavitation models at 𝜎 = 1.9. Figure 16 shows the time average void fraction field calculated with the incompressible (Incomp.) and compressible (Comp.) cavitation solvers and the experimental results. It can be seen that the simulations with two different cavitation solvers produce almost indistinguishable mean and RMS void fraction fields behind the wedge. The shape of mean cavity structures obtained numerically is comparable to the pattern observed in the experimental image. Note that there is a downstream offset in the mean and RMS void fraction fields between the numerical and the experimental results, which may be due to the low numerical resolution of the locally unstable flow structures within the detachment area. (a) (b) Figure 16. Comparisons of the experimentally (Exp.) [8] and numerically obtained time average (a) and (b) RMS void fraction fields using the incompressible (Incomp.) and the compressible (Comp.) cavitation models at 𝜎 = 1.9. Figure 16. Comparisons of the experimentally (Exp.) [ 8 ] and numerically obtained time average (a) and (b) RMS void fraction fields using the incompressible (Incomp.) and the compressible (Comp.) cavitation models at σ= 1.9. 6. Conclusions In this study, the numerical results obtained with both incompressible and compressible cavitation solvers were compared to assess the effects of fluid compressibility on the characteristics and dynamics of the cavitating flow behind a wedge. In relation to the experimental results obtained by Wu et al. [5], it was found: • Both cavitation solvers provide similar results to the experimental ones in terms of mean pressure and hydrodynamic forces. • Both cavitation solvers provide almost identical results of the dominant vortex shedding frequency and the instantaneous and mean void fraction fields. • The spectral content of the simulated hydrodynamic forces is similar with both solvers for low frequencies, but, for higher frequencies, the amplitudes are larger and the content is better resolved with the compressible solver. In conclusion, it was found that the compressibility effects on the cavitating vortex shedding behind the wedge can be simulated with the compressible solver and that they induce high-frequency phenomena, although the average quantities are not significantly affected in comparison with the incompressible solver results. Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/fluids9090218/s1, “Compressible cavitation solver.txt”: C code for the current compressible cavitation solver implemented in ANSYS Fluent. Author Contributions: Conceptualization, J.C. and L.G.; methodology, J.C.; software, J.C. and X.E.; validation, E.J. and X.E.; formal analysis, J.C. and L.G.; investigation, J.C. and L.G.; resources, E.J. and X.E.; data curation, E.J. and X.E.; writing—original draft preparation, J.C.; writing—review and editing, X.E.; visualization, E.J.; supervision, X.E.; funding acquisition, L.G. and X.E. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by Jiangsu Province Science Foundation for Youths, grant number BK20220538, and Jiangsu University, grant number 21JDG052. Fluids 2024,9, 218 19 of 20 Data Availability Statement: All the data generated or analyzed during this study are included in this article. Conflicts of Interest: The authors declare no conflicts of interest. References 1. Arndt, R.E.A. Cavitation in Fluid Machinery and Hydraulic Structures. Annu. Rev. Fluid Mech. 1981,13, 273–326. [CrossRef] 2. Escaler, X.; Egusquiza, E.; Farhat, M.; Avellan, F.; Coussirat, M. Detection of Cavitation in Hydraulic Turbines. Mech. Syst. Signal Process. 2006,20, 983–1007. [CrossRef] 3. Gu, Y.; Sun, H.; Wang, C.; Lu, R.; Liu, B.; Ge, J. Effect of Trimmed Rear Shroud on Performance and Axial Thrust of Multi-Stage Centrifugal Pump with Emphasis on Visualizing Flow Losses. J. Fluids Eng. 2024,146, 011204. [CrossRef] 4. Franc, J.P.; Michel, J.M. Fundamentals of Cavitation; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2006. 5. Young, J.O.; Holl, J.W. Effects of Cavitation on Periodic Wakes behind Symmetric Wedges. J. Basic Eng. 1966,88, 163–176. [CrossRef] 6. Belahadji, B.; Franc, J.P.; Michel, J.M. Cavitation in the Rotational Structures of a Turbulent Wake. J. Fluid Mech. 1995,287, 383–403. [CrossRef] 7. Ausoni, P.; Farhat, M.; Escaler, X.; Egusquiza, E.; Avellan, F. Cavitation Influence on von Kármán Vortex Shedding and Induced Hydrofoil Vibrations. J. Fluids Eng. 2007,129, 966–973. [CrossRef] 8. Wu, J.; Deijlen, L.; Bhatt, A.; Ganesh, H.; Ceccio, S.L. Cavitation Dynamics and Vortex Shedding in the Wake of a Bluff Body. J. Fluid Mech. 2021,917, A26. [CrossRef] 9. Nied´zwiedzka, A.; Schnerr, G.H.; Sobieski, W. Review of Numerical Models of Cavitating Flows with the Use of the Homogeneous Approach. Arch. Thermodyn. 2016,37, 71–88. [CrossRef] 10. Folden, T.S.; Aschmoneit, F.J. A Classification and Review of Cavitation Models with an Emphasis on Physical Aspects of Cavitation. Phys. Fluids 2023,35, 081301. [CrossRef] 11. Coutier-Delgosha, O.; Reboud, J.L.; Delannoy, Y. Numerical Simulation of the Unsteady Behaviour of Cavitating Flows. Int. J. Numer. Methods Fluids 2003,42, 527–548. [CrossRef] 12. Ji, B.; Luo, X.; Arndt, R.E.A.; Wu, Y. Numerical Simulation of Three Dimensional Cavitation Shedding Dynamics with Special Emphasis on Cavitation–Vortex Interaction. Ocean Eng. 2014,87, 64–77. [CrossRef] 13. Gnanaskandan, A.; Mahesh, K. Large Eddy Simulation of the Transition from Sheet to Cloud Cavitation over a Wedge. Int. J. Multiph. Flow 2016,83, 86–102. [CrossRef] 14. Callenaere, M.; Franc, J.P.; Michel, J.M.; Riondet, M. The Cavitation Instability Induced by the Development of a Re-Entrant Jet. J. Fluid Mech. 2001,444, 223–256. [CrossRef] 15. Kawanami, Y.; Kato, H.; Yamaguchi, H.; Tanimura, M.; Tagaya, Y. Mechanism and Control of Cloud Cavitation. J. Fluids Eng. 1997,119, 778–794. [CrossRef] 16. Laberteaux, K.R.; Ceccio, S.L. Partial Cavity Flows. Part 1. Cavities Forming on Models without Spanwise Variation. J. Fluid Mech. 2001,431, 1–41. [CrossRef] 17. Leroux, J.-B.; Astolfi, J.A.; Billard, J.Y. An Experimental Study of Unsteady Partial Cavitation. J. Fluids Eng. 2004,126, 94–101. [CrossRef] 18. Wu, J.; Ganesh, H.; Ceccio, S. Multimodal Partial Cavity Shedding on a Two-Dimensional Hydrofoil and Its Relation to the Presence of Bubbly Shocks. Exp. Fluids 2019,60, 66. [CrossRef] 19. Ganesh, H.; Mäkiharju, S.A.; Ceccio, S.L. Bubbly Shock Propagation as a Mechanism for Sheet-to-Cloud Transition of Partial Cavities. J. Fluid Mech. 2016,802, 37–78. [CrossRef] 20. Bhatt, A.; Ganesh, H.; Ceccio, S.L. Partial Cavity Shedding on a Hydrofoil Resulting from Re-Entrant Flow and Bubbly Shock Waves. J. Fluid Mech. 2023,957, A28. [CrossRef] 21. Gnanaskandan, A.; Mahesh, K. Numerical Investigation of Near-Wake Characteristics of Cavitating Flow over a Circular Cylinder. J. Fluid Mech. 2016,790, 453–491. [CrossRef] 22. Wang, C.; Wang, G.; Huang, B. Characteristics and Dynamics of Compressible Cavitating Flows with Special Emphasis on Compressibility Effects. Int. J. Multiph. Flow 2020,130, 103357. [CrossRef] 23. Vaca-Revelo, D.; Gnanaskandan, A. Numerical Assessment of the Condensation Shock Mechanism in Sheet to Cloud Cavitation Transition. Int. J. Multiph. Flow 2023,169, 104616. [CrossRef] 24. Katz, J. Cavitation Phenomena within Regions of Flow Separation. J. Fluid Mech. 1984,140, 397–436. [CrossRef] 25. O’Hern, T.J. An Experimental Investigation of Turbulent Shear Flow Cavitation. J. Fluid Mech. 1990,215, 365–391. [CrossRef] 26. Iyer, C.O.; Ceccio, S.L. The Influence of Developed Cavitation on the Flow of a Turbulent Shear Layer. Phys. Fluids 2002,14, 3414–3431. [CrossRef] 27. Choi, J.; Ceccio, S.L. Dynamics and Noise Emission of Vortex Cavitation Bubbles. J. Fluid Mech. 2007,575, 1–26. [CrossRef] 28. Agarwal, K.; Ram, O.; Lu, Y.; Katz, J. On the Pressure Field, Nuclei Dynamics and Their Relation to Cavitation Inception in a Turbulent Shear Layer. J. Fluid Mech. 2023,966, A31. [CrossRef] 29. Wang, Z.; Cheng, H.; Ji, B. Euler–Lagrange Study of Cavitating Turbulent Flow around a Hydrofoil. Phys. Fluids 2021,33, 112108. [CrossRef] Fluids 2024,9, 218 20 of 20 30. Wang, Z.; Cheng, H.; Ji, B.; Peng, X. Numerical Investigation of Inner Structure and Its Formation Mechanism of Cloud Cavitating Flow. Int. J. Multiph. Flow 2023,165, 104484. [CrossRef] 31. Brandao, F.L.; Bhatt, M.; Mahesh, K. Numerical Study of Cavitation Regimes in Flow over a Circular Cylinder. J. Fluid Mech. 2020, 885, A19. [CrossRef] 32. Park, S.; Seok, W.; Park, S.T.; Rhee, S.H.; Choe, Y.; Kim, C.; Kim, J.H.; Ahn, B.K. Compressibility Effects on Cavity Dynamics behind a Two-Dimensional Wedge. J. Mar. Sci. Eng. 2020,8, 39. [CrossRef] 33. Wang, Z.; Cheng, H.; Bensow, R.E.; Ji, B. Numerical Evaluation of the Bubble Dynamic Influence on the Characteristics of Multiscale Cavitating Flow in the Bluff Body Wake. Int. J. Multiph. Flow 2024,175, 104818. [CrossRef] 34. Zwart, P.; Gerber, A.G.; Belamri, T. A Two-Phase Flow Model for Predicting Cavitation Dynamics. In Proceedings of the Fifth International Conference on Multiphase Flow, Yokohama, Japan, 30 May–4 June 2004. 35. Egerer, C.P.; Schmidt, S.J.; Hickel, S.; Adams, N.A. Efficient Implicit LES Method for the Simulation of Turbulent Cavitating Flows. J. Comput. Phys. 2016,316, 453–469. [CrossRef] 36. Tong, S.-Y.; Zhang, S.; Wang, S.-P.; Li, S. Characteristics of the Bubble-Induced Pressure, Force, and Impulse on a Rigid Wall. Ocean Eng. 2022,255, 111484. [CrossRef] 37. Li, H.; Vasquez, S.A. Numerical Simulation of Steady and Unsteady Compressible Multiphase Flows. In Proceedings of the ASME 2012 International Mechanical Engineering Congress and Exposition, Houston, TX, USA, 9–15 November 2012; pp. 2239–2251. 38. Brunhart, M.; Soteriou, C.; Gavaises, M.; Karathanassis, I.; Koukouvinis, P.; Jahangir, S.; Poelma, C. Investigation of Cavitation and Vapor Shedding Mechanisms in a Venturi Nozzle. Phys. Fluids 2020,32, 083306. [CrossRef] 39. Mani, A. Analysis and Optimization of Numerical Sponge Layers as a Nonreflective Boundary Treatment. J. Comput. Phys. 2012, 231, 704–716. [CrossRef] 40. Schmidt, S.J. A Low Mach Number Consistent Compressible Approach for Simulation of Cavitating Flows. Ph.D. Thesis, Technical University of Munich, Munich, Germany, 2015. 41. Braza, M.; Chassaing, P.; Minh, H.H. Numerical Study and Physical Analysis of the Pressure and Velocity Fields in the near Wake of a Circular Cylinder. J. Fluid Mech. 1986,165, 79–130. [CrossRef] 42. Ding, H.; Shu, C.; Yeo, K.S.; Xu, D. Numerical Simulation of Flows around Two Circular Cylinders by Mesh-Free Least SquareBased Finite Difference Methods. Int. J. Numer. Methods Fluids 2007,53, 305–332. [CrossRef] 43. Seo, J.H.; Moon, Y.J.; Shin, B.R. Prediction of Cavitating Flow Noise by Direct Numerical Simulation. J. Comput. Phys. 2008,227, 6511–6531. [CrossRef] 44. Harichandan, A.B.; Roy, A. Numerical Investigation of Flow Past Single and Tandem Cylindrical Bodies in the Vicinity of a Plane Wall. J. Fluids Struct. 2012,33, 19–43. [CrossRef] 45. Qu, L.; Norberg, C.; Davidson, L.; Peng, S.-H.; Wang, F. Quantitative Numerical Analysis of Flow Past a Circular Cylinder at Reynolds Number between 50 and 200. J. Fluids Struct. 2013,39, 347–370. [CrossRef] 46. Kim, K.H.; Choi, J.I. Lock-in Regions of Laminar Flows over a Streamwise Oscillating Circular Cylinder. J. Fluid Mech. 2019,858, 315–351. [CrossRef] 47. Hong, S.; Son, G. Numerical Simulation of Cavitating Flows around an Oscillating Circular Cylinder. Ocean Eng. 2021,226, 108739. [CrossRef] 48. Menter, F.R. Best Practice: Scale-Resolving Simulations in ANSYS CFD; ANSYS Germany GmbH: Darmstadt, Germany, 2015. Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.