2.1. Cavity bubble inception

The fluid flowing through a contraction, result in the separation of the boundary layer and a substantial increase in velocity, accompanied by a turbulent field. The turbulence intensity significantly influences cavitation activity and its severity. Equation (2.1) presents a dimensionless number referred as cavitation number (CN), which quantifies the severity of cavitation within the flow domain

(2.1) \begin{align} CN=\frac{P_{d}-P_{v}} {P_{u}-P_{d}}, \end{align}

where P _d and P _v are the recovered downstream pressure and vapour pressure of the liquid at constant temperature, respectively, P _u is the upstream pressure in fluid domain. The CN provides the intensity of cavitation phenomenon in the fluid domain due to a contraction in flow geometry which leads to an increase in the velocity of the flow lamina. The higher pressure gradient across the flow domain results in a lower value of CN, which means more inception of the cavity bubbles and higher chances of material removal (Mehdi Hadi, Reference Hadi2013). Thus, the cavity bubble initiation occurs at CN <1. The governing equation and mathematical formulation to simulate the fluid flow and cavitation phenomenon are described further.

2.2. Turbulence model

In this section the governing models for a homogeneous multiphase (water and vapour) cavitating waterjet through an orifice are discussed. A realisable k − ε turbulence model with standard wall function is used. This model, based on Reynolds-averaged NavierStokes (RANS) equations, effectively addresses mathematical instability inherent in highly turbulent flows. In the present study, the homogeneous mixture of water and vapour is considered as a single fluid. Consequently, the homogeneous mixture fluid is appropriately coupled to mass transfer equations. The primary governing equations, along with the turbulence model for CFD simulation, are elaborated in the subsequent sub-section. The continuity equation of homogeneous mixture is expressed in equation (2.2)

(2.2) \begin{align} \frac{\partial \rho _{m}} {\partial t}+\nabla \rho _{m} \vec {v}_{m}=0, \end{align}

where ρ _m and $ \vec {v}_{m} $ are the density and velocity of the mixture phase, respectively. Assuming both phases flow at the same velocity, the momentum equation for the homogeneous mixture fluid is given in equation (2.3)

(2.3) \begin{align}\frac {\partial \left(\rho _{m}\vec {v}_{m}\right)}{\partial t}+\nabla \cdot \left(\rho _{m}\vec {v}_{m}\vec {v}_{m}\right)=-\nabla p+\nabla \cdot \left[\mu _{m}(\nabla \vec {v}_{m}+\nabla {\vec {v}_{m}}^{T})\right]+\rho _{m}\overset{\tiny\rightharpoonup}{g}+\overset{\tiny\rightharpoonup}{F}, \end{align}

where μ _m is the viscosity of the mixture fluid, $ \rho_{m} \overset{\tiny\rightharpoonup}{g} $ is the gravitational body force and the term $ \overset{\tiny\rightharpoonup}{F} $ accounts for the added body forces applied to the mixture fluid volume (this force may arise from the interaction with a secondary phase). Therefore, the Boussineq hypothesis (Guide, 2013) is used to close the RANS equation. In addition, k is the turbulence kinetic energy and ε is the turbulence dissipation rate which are modelled using the transport equations (2.4)–(2.5)

(2.4) \begin{align} & \frac {\partial}{\partial t}\left(\rho k\right)+ \frac{\partial}{\partial x_{i}}\left(\rho kv_{i}\right)= \frac{\partial}{\partial x_{i}}\left(\Gamma _{k} \frac{\partial k }{\partial x_{i}}\right)+G_{k}-G_{b}-\rho \varepsilon -Y_{k}+S_{k}, \end{align}

(2.5) \begin{align} & \frac{\partial}{\partial t}\left(\rho \varepsilon \right) + \frac{\partial}{\partial x_{i}}\left(\rho \varepsilon v_{i}\right) = \frac{\partial }{ \partial x_{i}}\left(\Gamma _{\varepsilon } \frac{\partial \varepsilon }{\partial x_{i}}\right)+\rho C_{1}S\varepsilon -\rho C_{2} \frac{\varepsilon ^{2}}{k+\sqrt{\nu \varepsilon }}+C_{1\varepsilon } \frac{\varepsilon}{k}C_{3\varepsilon }G_{b}+S_{\varepsilon }, \\[10pt] \nonumber \end{align}

where $C_{1}=\mathit{\max } \left[0.43, \frac{\eta }{\eta +5}\right],\eta =S \frac{k }{\varepsilon },S=\sqrt{2S_{ij}S_{ij}}$

In equations (2.4)-(2.5), G _k is the initiation of k due to the mean velocity gradient. G _b is used for the generation of k due to the buoyancy effect. However, in the present study the buoyancy effect on ε is neglected, by placing the value of G _b to zero in equation (2.5). Also, Γ _k and Γ _ε represent the effective diffusivity of k and ε, respectively, Y _k is the dissipation of k due to turbulence, C ₂ and C _1ε are constant values of 1.9 and 1.44, respectively (Guide, 2013), and S _k and S _ε are user defined source terms. The cavitation model is discussed in the following sub-section.

2.3. Cavitation model

The cavitation model proposed by Schnerr and Sauer (Schnerr and Sauer, Reference Schnerr and Sauer2001) has been used in the present study. This model helps in determining flow properties of the mixture fluid as a linear function of vapour volume fraction. Therefore, the density of the mixture fluid is shown in equation (2.6)

(2.6) \begin{align} \rho _{m}=\alpha \rho _{v}+\left(1-\alpha \right)\rho _{l}. \end{align}

The liquid–vapour mass transfer occurs in the limited zone for incompressible homogeneous mixture flow, when the static pressure is lower than or equal to the vapour pressure at constant temperature of the working fluid. The vapour transport equation that governs the mass transfer between liquid and vapour is given in equation (2.7)

(2.7) \begin{align} \frac {\partial}{\partial t}\left(\alpha \rho _{v}\right)+\nabla \cdot \left(\alpha \rho _{v}\vec {V}\right)=R, \end{align}

where α stands for vapour fraction, ρ _v stands for density of the vapour phase, $ \vec{V} $ is the velocity. Equation (2.8) shows the net mass source term R from the liquid to vapour

(2.8) \begin{align} R= \frac{\rho _{v}\rho _{l}}{\rho }\cdot \frac{d\alpha}{dt}. \end{align}

The cavitation models considered in the present research are derived based on various assumptions such as the mixture is supposed to be incompressible and homogeneous; the nuclei that already exist in the flow regime are the source for the cavity bubble, which can grow and implode depending on the surrounding conditions (temperature and pressure); the slip velocity between the phases is neglected; positive mass transfer from liquid to vapour is considered; and the bubble growth is represented by a simplified Rayleigh–Plesset equation.

Nonetheless, the carrier fluid contains a large number of nuclei that are needed for the cavity bubble to form. Therefore, the simplified Rayleigh–Plesset equation, (2.9), provides estimation of bubble growth and collapse in a flowing liquid (Plesset and Winet, Reference Plesset and Winet1974)

(2.9) \begin{align} \frac{dR_{B}}{dt}=\sqrt{\left( \frac{2}{3} \frac{P\left(R_{B}\right)-P_{\infty } }{\rho _{l}}\right)}. \end{align}

Initially, we starts with the bubble nucleus of radius R _B and density per unit volume of liquid n _o . The relationship between the vapour volume fractions to the number of bubble nuclei for a unit volume of mixture fluid is shown in equation (2.10)

(2.10) \begin{align} \alpha _{vo}=\frac{n_{o}\cdot \frac{4}{3}\pi R_{B}^{3}}{1+n_{o}\frac{4}{3}\pi R_{B}^{3}}. \end{align}

The relationship between bubble radius, volume fraction of phase and bubble number density is expressed in equation (2.11)

(2.11) \begin{align} R_{B}=\left(\frac{\alpha _{v}}{\alpha _{l}}\frac{3}{4\pi n_{o}}\right)^{1/3}. \end{align}

Schnerr and Sauer (Schnerr and Sauer, Reference Schnerr and Sauer2001) proposed the modified mass transfer model by considering equations (2.6–2.11). The modified mass transfer model is defined as in equation (2.12)

(2.12) \begin{align} \dot{m}_{e}^{\prime\prime\prime}= \frac{\rho _{l}\rho _{v}}{\rho }\alpha _{l}\alpha _{v}\frac{3}{R_{B}}\sqrt{\frac{2}{3} \frac{P\left(R_{B}\right)-P_{\infty }}{\rho _{l}}}. \end{align}

Now, n _o and R _B are input parameters for the modified mass transfer model and the standard values are 10⁸ and 10^-5 m, respectively. Thereafter, the next section presents the Lagrangian tracking of a bubble in the cavitation zone.

2.4. Lagrangian tracking of bubbles

In Discrete Phase Model (DPM) in Fluent also termed as Discrete Bubble Model (DBM), Lagrangian approach is used to track the discrete phase (cavity bubbles) in fluid domain, whereas the continuous phase (water) is solved using the transport equation for mass, momentum, and turbulent models. The discrete phase is considered as a point. Each point has its own position, velocity, radius, bubble pressure, etc. Considering these properties bubble-bubble interaction models (coalescence, breakage), bubble–wall interactions and gas diffusion at the bubble interface can be modelled. As the objective of the present research is to observe the phenomenon of CWJM, one-way coupling is used. One-way coupling is applied for interactions between the phases (continuous and discrete). Thus, the mass and momentum equations for the continuous phase are solved in single-phase form, with the resulting velocity and pressure of the continuous phase used to determine the trajectory of the discrete phase by applying additional forces. The mathematical formulation of the models are discussed in details in the authors’ previous publication (Kumar et al., Reference Kumar, Bera and Rout2024).

2.5. Bubble diameters and statistics

The Lagrangian tracking of cavity bubbles in cavitating flow offers the advantage of individually tracking each bubble, allowing for more precise applications for bubble dynamics. The vapour bubbles of random diameters are injected from the inlet surface. The vapour bubbles are initialised with the fluid velocity in Eulerian reference. Both the phases are in physical, chemical and thermal equilibrium with activated breakage and coalescence models. There are many experimental methods to measure the size of bubbles in the flow domain such as the light transmission method, Mastersizer (Couto et al., Reference Couto, Nunes, Neumann and França2009), focused beam reflectance measurement (Kail et al., Reference Kail, Marquardt and Briesen2009) and acoustic-based spectrometry (Wu and Chahine, Reference Wu and Chahine2010). But, it was concluded experimentally, due to low conductance of the gas phases, noise and the interfering signal emitted from the bubbles, these methods failed to predict accurate results.

Therefore, the Sauter mean diameter (SMD) calculates the average bubble size in the flow domain (Kowalczuk and Drzymala, Reference Kowalczuk and Drzymala2016). This parameter, introduced by Sauter in 1926 in his study of fuel oil droplets, is referred as surface–volume mean. The calculation of SMD of different sizes of bubbles considering both volume and surface area of the bubbles is presented in equation (2.13)

(2.13) \begin{align} SMD\left(d_{32}\right)= \frac{\sum _{i=1}^{n}n_{i}d_{i}^{3}}{\sum _{i=1}^{n}n_{i}d_{i}^{2}}, \end{align}

where n _i and d _i are the number and diameter of the individual cavity bubbles, respectively. Currently, SMD is widely recognised for determining the average size of gas bubbles and liquid droplets.

2.6. Impact energy modelling

The surface damage caused by hydrodynamic cavitation on materials primarily occurs in the form of micro-pits, resulting from micro-plastic deformation due to the impact of microjets. Franc (Franc, Reference Franc2009) performed a phenomenological study of the cavitation erosion process on a ductile material, they shed light on two characteristics scales such as the “time scale” which is the time required for impacts to completely cover the specimen surface and the “length Scale”, this is the critical length for a specimen that governs the erosion behaviour. They also proposed that the incubation time is proportional to the covering time and influenced by the aggressiveness of the flow. The erosion rate under steady-state conditions is determined by the ratio of the hardened layer thickness to the covering time. In the incubation period, the impact pressure resulting from the collapse of a cavitation bubble is complex in both spatial and temporal dimensions, exhibiting unsteady characteristics such as the formation of microjets and shock waves (Brennen, Reference Brennen1995). Hence, a more effective approach for the inverse finite element method (FEM) would be to utilise the impact load generated by the collapse of cavity bubbles. However, for this FEM a Fluid-Structure Interaction (FSI) method must be considered and the authors in Hsiao et al. (Reference Hsiao, Jayaprakash, Kapahi, Choi and Chahine2014) and Roy et al. (Reference Roy, Franc, Pellone and Fivel2015) have pointed out the limitations of such simulations. The only possible way is to use a representative pressure field and analyse how the materials respond to the impact load.

Therefore, an analytical procedure is considered where the shape of the pit is in round shape like a bowl, as shown in Figure 3. The assumed size of the pit from implosion of a cavity bubble is completely characterised by its depth (h _p ) and diameter (d _p ). It is common to assume a pit to be spherical in shape, if the mean aspect ratio of its diameter to depth is large (Reference Carnelli, Karimi and Franc2012a, Carnelli et al., Reference Carnelli, Karimi and Franc2012b; Tzanakis et al., Reference Tzanakis, Eskin, Georgoulas and Fytanidis2014). The assessment of pits can be carried out by analysing the geometrical shape of the deformed part and the plastic strain induced in the region. An analogy is developed between the material deformation caused by microjet liquid impact and spherical indentation. Considering the scope of the current study, the fundamental equations of the spherical indentation tests are considered.

Figure 3. Schematic representation of the spherical geometry used to model cavitation pit, including the extension of the plastically deformed area behind the pit.

Initially, Tabor (Tabor, Reference Tabor1956) has presented the indentation response on the metallic workpiece. Using the definition of true indentation stress (σ) on a metal workpiece, equation (2.14) can be written as

(2.14) \begin{align} \sigma = \frac{P_{m}}{\psi }= \frac{1}{\psi }\frac{L}{A}, \end{align}

where P _m and L are the mean contact pressure and indentation load on metal workpiece, respectively, over the projected area A and ψ is a constraint parameter introduced by Francis (Francis, Reference Francis1976). It was noted that ψ depends on the indentation regime occurring in the core beneath the indenter as expressed in equation (2.15)

(2.15) \begin{align} \psi =\left\{\begin{array}{c@{\quad}c} 1.1 & \phi \lt 1.1\\ 1.1+0.5log\phi & 1.1\lt \phi \lt 27.3\\ 2.87 & \phi \gt 27.3 \end{array}\right\}, \end{align}

where $ \phi $ is a dimensionless parameter associated with the elastic strain at the onset of deformation in the metallic workpiece. If the workpiece is in a purely elastic regime, the constraint factor is 1.1. As the workpiece transitions to an elastic-plastic regime, the constraint factor increases, reaching a maximum value of 2.87 when the regime becomes fully plastic. The representative plastic strain ε and the indentation parameters as proposed by Tabor (Tabor, Reference Tabor1956), are presented in equation (2.16)

(2.16) \begin{align} \epsilon =0.2\frac{d}{D}, \end{align}

where $ \epsilon $ represents the uniaxial strain at the contact edge of the indenter, and d is the diameter of the spherical cap impression created by the spherical indenter having diameter D. The nonlinear relationship between stress and strain, as well as the work hardening occurring on the surface of the workpiece, is modelled using Hollomon’s power law and is represented in equation (2.17)

(2.17) \begin{align} \sigma =\left\{\begin{array}{cc} \begin{array}{c} \varepsilon E\\ \sigma _{y}\left( \frac{\varepsilon}{\varepsilon _{y}}\right)^{n} \end{array} & \begin{array}{c} \varepsilon \leq \varepsilon _{y}\\ \varepsilon \gt \varepsilon _{y} \end{array} \end{array}\right\}, \end{align}

where σ _y represents the yield stress at reference stress at 0.002 of plastic strain for the specific workpiece material, ε represents the total strain (elastic and plastic), ε _y is the strain at yield point. In the present work the strain value is taken as 0.002, and n is the strain hardening exponent that describes the hardening behaviour of the workpiece sample, which depends on the materials of the workpiece (Bonnavand et al., Reference Bonnavand, Bramley and Mynors2001).

The above fundamental relationship is used to establish the relationship of the pit stress analysis from the microjet impact during implosion of the cavitation bubbles. Therefore, a reverse engineering approach is used to establish the relation between the material deformation with the impact of the microjet and indentation (Carnelli et al., Reference Carnelli, Karimi and Franc2012a). The deformation of materials due to the impingement on implosion of an individual sphere-shaped cavity bubble at the surface of the workpiece is known as the cavitation pit. Assuming that a cavitation pit could be like a bowl geometry, and completely described by its depth (h _p ) and diameter (d _p ) for modelling the hydrodynamic impact on the workpiece. The magnitude of pit depth and diameters are used to determine strain at the boundary of the cavitation pit. Considering equation (2.16), the equivalent uniaxial strain ε _p at the contact boundary of the spherical cup-like pit can be expressed as in equation (2.18)

(2.18) \begin{align} \varepsilon _{p}=0.2\frac{d}{D}= \frac{\frac{d_{p}}{2}}{\sqrt{\left(\frac{d_{p}}{2}\right)^{2}+\left(R-h_{p}\right)^{2}}}=0.2\frac{h_{p}d_{p}}{\left(\frac{d_{p}}{2}\right)^{2}+h_{p}^{2}}, \end{align}

where d _p and h _p represent the deformed pit diameter and its depth, respectively, with a cavity bubble of radius R. The induced stressσ _p within the cavitation pit can be determined from the stress–strain relationship outlined in equation (2.17), in conjunction with equation (2.18). Therefore, σ _p is expressed in equation (2.19)

(2.19) \begin{align} \sigma _{p}=\sigma _{y}\left(\frac{\varepsilon _{p}}{\varepsilon _{y}}\right)^{n}. \end{align}

After determining the stress induced on the cavitation pit due to the impact of the microjet on the surface of the workpiece, the impact load L _p as in equation (2.20) can be derived from the relationship shown in equation (2.14)

(2.20) \begin{align} L_{p}=\sigma _{p}\psi A_{p}. \end{align}

After calculating the impact load, the final step is to determine the hydrodynamic impact pressure (σ _H ) produced during the implosion of the cavity bubble. This is derived by considering the geometrical characteristics of the needle-like pits and the induced stress. For calculation of hydrodynamic impact stress (σ _H ), it is assumed that the load L _p measured from the pit stress is equivalent to the hydrodynamic load L _H released during the implosion of the cavity bubble. It is basically the impact load generated on surface of the workpiece due to the microjet and the expression is shown in equation (2.21)

(2.21) \begin{align} L_{H}=\int _{0}^{R}\mathit{\sigma }SRdR=\sigma _{H}\pi \frac{D_{H}^{2}}{4}, \end{align}

where S represents the perimeter of the impacted area and D _H is the effective hydrodynamic distance, defined as the diameter of the microjet in this study. The diameter of the microjet is equal to 1/10th of the resonant bubble diameter, as explained in previous studies (Tezel and Mitragotri, Reference Tezel and Mitragotri2003). Furthermore, the theory of bubble collapse proposed by Rayleigh (Rayleigh, Reference Rayleigh1917) is considered where the maximum energy (E _b ) released due to implosion of the cavity bubble is expressed as in equation (2.22)

(2.22) \begin{align} E_{b}=\frac{4}{3}\pi \left(P_{\textit{static}}-P_{\textit{vapour}}\right)R_{\textit{resonant}}^{3}, \end{align}

where P _static and P _vapour represent the static pressure and vapour pressure inside the bubble (3500 Pa at 27 °C), respectively, and R _resonant is the resonant radius of the bubble. The energy equation states that the energy of a cavity bubble is determined by its maximum volume and the associated pressure gradient. The energy generated by the bubble is then transformed into significant kinetic energy, propelling the microjet toward the surface of the workpiece. The hydrodynamic impact of microjet causes a high pressure and its magnitude can be obtained from the water hammer of equation (2.23) (Tijsseling and Anderson, 2004)

(2.23) \begin{align} \sigma _{H}=\rho _{l}C_{l}\left(\frac{\rho _{s}C_{s}}{\rho _{l}C_{l}+\rho _{s}C_{s}}\right)V_{jet}, \end{align}

where ρ _l and C _l represent the density and speed of sound in the liquid medium. Similarly, ρ _s and C _s are the relating parameters for the workpiece.

2.7. Numerical simulation set-up

The model was implemented within the commercial CFD software package ANSYS Fluent, version 18.1. The complete simulation process is illustrated in Figure 4, which begins with the creation of a geometric model, followed by mesh generation and boundary definition. Subsequently, the solution convergence is verified, and upon achieving convergence, the results are visualised. This is followed by selection of the output parameters for the grid independence test to ensure the desired accuracy while minimising computation time.

Figure 4. Multiphase CFD simulation approach.

The numerical method for simulating cavitation phenomenon in the flow domain starts with initialising the solution using specific boundary conditions and fluid properties. This is followed by solving the RANS equation, turbulence model and cavitation model. The mass transfer equation is solved to capture the interaction between the bulk liquid and discrete phases. Once convergence criteria are met, these equations are iteratively solved for the n time step.