15. Homogeneous Cavitation Modeling

Cavitation is the formation of vapor when the local pressure becomes low enough for liquid to evaporate. Vapor can condense when the pressure rises. The large difference between liquid and vapor densities couples this mass transfer strongly to the pressure and velocity fields. Evaporation also absorbs heat, changing temperature and therefore saturation pressure.

This chapter develops the formulation used by cavitation_nist: a compressible liquid–vapor mixture with optional noncondensable gas (NCG), finite-rate phase change, and a temperature equation. It includes the optional transported bubble population and turbulence corrections. General flow discretization and pressure–velocity coupling follow the preceding chapters.

15.1. Physical Model and Assumptions

Liquid, vapor, and NCG share one velocity \(\vec u\), pressure \(p\), and temperature \(T\). Subscripts \(l\), \(v\), and \(g\) identify these constituents. Liquid and vapor are phases of the same condensable substance. NCG is transported with the mixture but does not participate in evaporation or condensation.

The homogeneous approximation neglects phase slip and separate phase temperatures. Bubble interfaces are not resolved. Their growth or collapse enters through a volumetric mass-transfer model, which determines how fast composition changes. Common pressure and temperature do not require the vapor content to equal its equilibrium value.

Thermal coupling is retained because evaporation cools the mixture and can reduce the saturation pressure that drives further evaporation. This effect can be appreciable in cryogenic liquids. The model represents it through the mixture temperature and temperature-dependent properties; it does not resolve the thermal layer around each bubble.

15.2. Phase Fractions and Thermodynamic Closure

15.2.1. Mass, Volume, and Density

Let \(Y_i\) be mass fraction and \(\alpha_i\) volume fraction. Admissible composition satisfies

(15.2.1.1)\[Y_l=1-Y_v-Y_g \qquad Y_i\geq0 \qquad \sum_iY_i=1\]

An intrinsic phase density \(\rho_i(p,T)\) is mass per volume of that phase; the partial density \(\rho Y_i\) is mass per mixture volume. Adding the phase volumes gives

(15.2.1.2)\[\frac{1}{\rho}=\sum_{i=l,v,g}\frac{Y_i}{\rho_i(p,T)} \qquad \alpha_i=\frac{\rho Y_i}{\rho_i(p,T)} \qquad \sum_i\alpha_i=1\]

The general additive-volume construction is given in Additive-Volume Real-Fluid EOS Closure. Here the liquid and vapor densities come from separate real-fluid property tables at the current \(p,T\). For the calorically perfect NCG,

(15.2.1.3)\[\rho_g=\frac{p}{R_gT} \qquad e_g=(c_{p,g}-R_g)T \qquad h_g=c_{p,g}T\]

Here \(R_g\) and \(c_{p,g}\) are constant gas properties. The energy datum in Eq. (15.2.1.3) is consistent with the gas closure; the liquid and vapor retain their tabulated energy datums.

15.2.2. Phase Properties and Saturation Properties

The saturation table supplies \(p_{\mathrm{sat}}(T)\), saturated liquid and vapor properties, and surface tension \(\sigma(T)\). Use \(\rho_i^{\mathrm{sat}}(T)\) for saturation densities. They differ from \(\rho_i(p,T)\) away from saturation. A superscript \(\mathrm{sat}\) will therefore be retained in the source formulas that use them.

Property interpolation is linear in temperature at each pressure level, then linear between pressure levels. Queries outside the available range use endpoint properties. This numerical extension does not establish physical validity beyond the table. Saturation data are also interpolated in temperature, and the phase-change sources are inactive at and above the upper saturation temperature. Below the lower saturation temperature, the source properties use their lower endpoint values.

At a boundary with specified pressure, density, and composition, temperature is recovered from the complete mixture EOS as described in Additive-Volume Real-Fluid EOS Closure. Incoming flow uses its prescribed composition; outgoing flow uses extrapolated composition. The NCG term must be retained because even a small gas mass can occupy appreciable volume.

15.2.3. Mixture Properties

The caloric closure assumes mass-weighted specific energy, with no excess energy of mixing. Together with volume additivity, this gives mass-weighted enthalpy and heat capacity:

(15.2.3.1)\[e=\sum_iY_ie_i \qquad h=\sum_iY_ih_i \qquad c_p=\sum_iY_ic_{p,i}\]

The thermal expansion coefficient of a phase is \(\beta_{T,i}=-(1/\rho_i)(\partial\rho_i/\partial T)_p\). Differentiating Eq. (15.2.1.2) at fixed pressure and composition gives

(15.2.3.2)\[\beta_{T,m} =-\frac{1}{\rho}\left(\frac{\partial\rho}{\partial T}\right)_{p,Y} =\sum_i\alpha_i\beta_{T,i}\]

The temperature equation uses the dimensionless factor \(\beta_{T,m}T\). For ideal-gas NCG, \(\beta_{T,g}T=1\). Molecular viscosity and thermal conductivity use volume weights:

(15.2.3.3)\[\mu=\sum_i\alpha_i\mu_i \qquad \lambda=\sum_i\alpha_i\lambda_i\]

The adopted acoustic mixture rule uses the pure-phase sound speeds \(a_i\):

(15.2.3.4)\[\frac{1}{\rho^2a^2}=\sum_i\frac{Y_i}{\rho_i^2a_i^2} \qquad a_g^2=\frac{c_{p,g}}{c_{p,g}-R_g}R_gT\]

All phase properties in Eqs. (15.2.3.1)–(15.2.3.4) are evaluated at the common current state. Equation (15.2.3.4) approximates acoustic propagation at fixed composition. It excludes the additional response from phase change during a pressure disturbance.

15.3. Conservation Laws and Thermal Coupling

15.3.1. Mixture and Phase Mass Balances

Let \(\dot m_{lv}\) be mass transferred from liquid to vapor per mixture volume per time, in \(\mathrm{kg\,m^{-3}\,s^{-1}}\). It is positive for evaporation and negative for condensation. With no external mass injection,

(15.3.1.1)\[\begin{split}\begin{aligned} \frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\vec u)&=0\\ \frac{\partial(\rho Y_v)}{\partial t} +\nabla\cdot(\rho Y_v\vec u)&=\dot m_{lv}\\ \frac{\partial(\rho Y_g)}{\partial t} +\nabla\cdot(\rho Y_g\vec u)&=0 \end{aligned}\end{split}\]

Subtracting the vapor and gas equations from continuity gives the liquid source \(-\dot m_{lv}\). Thus phase change redistributes mass without creating a mixture mass source. These equations contain no separate phase diffusion or dissolved-gas exchange model.

The momentum equation is the single-mixture equation already developed in the guide, with the density and viscosity defined above. Interphase momentum exchange does not require a separate equation because velocity is shared. Bubble surface tension can affect a growth law without adding a resolved-interface capillary force to this mixture momentum equation.

15.3.2. From Enthalpy to Temperature

Write the static-enthalpy balance as

(15.3.2.1)\[\rho\frac{Dh}{Dt}=\frac{Dp}{Dt}+\mathcal H \qquad \mathcal H=-\nabla\cdot\vec q +\Matrix{\tau}:\nabla\vec u+Q_h\]

Here \(D/Dt=\partial/\partial t+\vec u\cdot\nabla\), \(\vec q\) is the modeled heat flux, \(\Matrix{\tau}\) is viscous stress, and \(Q_h\) denotes other prescribed volumetric heating. Turbulent heat transport and stress closures follow RANS-Based Turbulence Models.

For each phase, the thermodynamic enthalpy differential is

(15.3.2.2)\[dh_i=c_{p,i}\,dT+\frac{1-\beta_{T,i}T}{\rho_i}\,dp\]

Differentiate \(h=\sum_iY_ih_i\) and eliminate \(dY_l\) using \(dY_l=-dY_v-dY_g\). Equations (15.2.1.2) and (15.2.3.2) then give

(15.3.2.3)\[dh=c_p\,dT+\frac{1-\beta_{T,m}T}{\rho}\,dp +(h_v-h_l)\,dY_v+(h_g-h_l)\,dY_g\]

Continuity and the phase balances imply \(\rho DY_v/Dt=\dot m_{lv}\) and \(DY_g/Dt=0\). Substitution into Eq. (15.3.2.1) yields

(15.3.2.4)\[\rho c_p\frac{DT}{Dt} =\mathcal H+\beta_{T,m}T\frac{Dp}{Dt} -(h_v-h_l)\dot m_{lv}\]

The cavitation thermal closure approximates the phase enthalpy difference by the saturation latent heat at the mixture temperature:

(15.3.2.5)\[L(T)=h_v^{\mathrm{sat}}(T)-h_l^{\mathrm{sat}}(T)\]

With this approximation, the temperature equation is

(15.3.2.6)\[\boxed{ \rho c_p\frac{DT}{Dt} =\mathcal H+\beta_{T,m}T\frac{Dp}{Dt}-L(T)\dot m_{lv} }\]

Evaporation removes sensible heat when \(L>0\); condensation releases it. The pressure-work multiplier is generally not one for a real-fluid mixture. Equation (15.3.2.4) follows from differentiable, consistent phase properties; Eq. (15.3.2.6) additionally makes the saturation latent-heat approximation. Separately interpolated tables need not satisfy every thermodynamic derivative identity exactly.

15.4. Evaporation and Condensation Models

Use \(p_s\) for the pressure threshold in the source law. Normally \(p_s=p_{\mathrm{sat}}(T)\); a prescribed threshold or the turbulent correction described below can replace it. Use \([x]_+=\max(x,0)\). The net rate is evaporation minus condensation. These source laws apply below the upper saturation temperature.

15.4.1. Merkle Model

The Merkle model scales the pressure departure by a reference dynamic pressure and the mass-transfer time by a reference flow time [MFB98]. Let \(U_r\), \(\ell_r\), and \(\rho_r\) be reference velocity, length, and density, and define

(15.4.1.1)\[t_r=\frac{\ell_r}{U_r} \qquad q_r=\frac12\rho_rU_r^2 \qquad k_+=\frac{[p_s-p]_+}{\tau_+t_rq_r} \qquad k_-=\frac{[p-p_s]_+}{\tau_-t_rq_r}\]

The dimensionless positive parameters \(\tau_+\) and \(\tau_-\) set the evaporation and condensation scales. Both \(k_+\) and \(k_-\) have units \(\mathrm{s^{-1}}\). The source is

(15.4.1.2)\[\dot m_{lv}=\rho\left(k_+Y_l-k_-Y_v\right)\]

Available liquid limits evaporation and available vapor limits condensation. The coefficients are empirical; the reference scales and their calibration are part of the model.

15.4.2. Bubble-Growth Approximation

For a spherical bubble in an incompressible liquid, a Rayleigh–Plesset description is [Bre95]

(15.4.2.1)\[R_b\ddot R_b+\frac32\dot R_b^2 =\frac{p_s+p_b-p}{\rho_l} -\frac{4\mu_l\dot R_b}{\rho_lR_b} -\frac{2\sigma}{\rho_lR_b}\]

Here \(p_b\) is the gas partial pressure inside the bubble, and \(p\) is the surrounding liquid pressure. The radius derivatives follow the bubble. Neglecting radial acceleration, viscosity, surface tension, and gas pressure gives \(\dot R_b^2=2(p_s-p)/(3\rho_l)\) for evaporation. Using saturation liquid density and prescribing a corresponding collapse speed gives the positive branch speeds

(15.4.2.2)\[u_+=\sqrt{\frac{2[p_s-p]_+}{3\rho_l^{\mathrm{sat}}}} \qquad u_-=\sqrt{\frac{2[p-p_s]_+}{3\rho_l^{\mathrm{sat}}}}\]

The signed radial speed is \(\dot R_b=u_+-u_-\). The collapse branch is a prescribed continuation: the acceleration-free quadratic alone has no real root for \(p>p_s\).

For fixed vapor density, growth of one bubble changes its vapor mass at the rate \(4\pi R_b^2\rho_v\dot R_b\). Multiplication by the number of bubbles per mixture volume supplies the starting point for the bubble-based mass-transfer closures.

15.4.3. Zwart Model

The Zwart model uses a prescribed bubble radius \(R_0\), a nuclei volume parameter \(\alpha_0\), and empirical coefficients \(F_+\) and \(F_-\) [ZGB04]. Define the source-model vapor and residual liquid volumes by

(15.4.3.1)\[\widehat\alpha_v=\frac{\rho Y_v}{\rho_v^{\mathrm{sat}}} \qquad \widehat\alpha_l=1-\widehat\alpha_v-\frac{\rho Y_g}{\rho_g}\]

The nuclei concentration is represented by \(\alpha_0\widehat\alpha_l\) during evaporation. During condensation, the existing vapor volume is used. The resulting source is

(15.4.3.2)\[\dot m_{lv} =\frac{3\rho_v^{\mathrm{sat}}}{R_0} \left(F_+\alpha_0\widehat\alpha_l u_+ -F_-\widehat\alpha_v u_-\right)\]

The hat distinguishes these source estimates from the EOS volume fractions \(\alpha_i\). They coincide with the corresponding physical fractions when the relevant phase densities agree. Away from saturation, or with lagged mixture density, \(\widehat\alpha_l\) is not guaranteed to be a nonnegative physical liquid fraction. This is a limitation of this source closure.

15.4.4. Fixed-Count Schnerr–Sauer Model

In the two-phase fixed-count construction, let \(n_0\) be the number of bubbles per liquid volume. For equal-radius bubbles, \(V_v=n_0V_l(4\pi R_b^3/3)\), so

(15.4.4.1)\[\frac{\alpha_v}{1-\alpha_v}=\frac{4\pi}{3}n_0R_b^3 \qquad R_b=\left(\frac{3\alpha_v}{4\pi n_0(1-\alpha_v)}\right)^{1/3}\]

For constant phase densities, \(\rho=\rho_l(1-\alpha_v)+\rho_v\alpha_v\). Differentiating \(Y_v=\rho_v\alpha_v/\rho\) and the radius relation at fixed \(n_0\) gives

(15.4.4.2)\[\frac{dY_v}{d\alpha_v}=\frac{\rho_l\rho_v}{\rho^2} \qquad \frac{D\alpha_v}{Dt} =3\alpha_v(1-\alpha_v)\frac{\dot R_b}{R_b}\]

Substitute these relations into the vapor balance to obtain the two-phase Schnerr–Sauer source [SS01]:

(15.4.4.3)\[\dot m_{lv} =\rho\frac{dY_v}{d\alpha_v}\frac{D\alpha_v}{Dt} =\frac{3\rho_l\rho_v}{\rho}\alpha_v(1-\alpha_v)\frac{\dot R_b}{R_b} =3\rho Y_vY_l\frac{\dot R_b}{R_b}\]

The fixed-count closure uses saturation densities and a seed radius \(R_0\). Its inception seed and effective vapor fraction are

(15.4.4.4)\[z_0=\frac{4\pi}{3}n_0R_0^3 \qquad \alpha_0=\frac{z_0}{1+z_0} \qquad Y_0=\frac{\rho_v^{\mathrm{sat}}\alpha_0}{\rho} \qquad Y_e=\max(Y_v,Y_0)\]

Define \(\widehat\alpha_e=\rho Y_e/\rho_v^{\mathrm{sat}}\). The radius and source used with this construction are

(15.4.4.5)\[\begin{split}\begin{aligned} R_b&=\left(\frac{3\widehat\alpha_e} {4\pi n_0(1-\widehat\alpha_e)}\right)^{1/3}\\ \dot m_{lv}&=\frac{3\rho_l^{\mathrm{sat}}}{R_b} \widehat\alpha_l\left(Y_eu_+-Y_vu_-\right) \end{aligned}\end{split}\]

The radius relation requires \(\widehat\alpha_e<1\); a small positive denominator floor regularizes its numerical evaluation. The seed permits inception from zero vapor without adding transported phase mass. NCG enters the liquid-availability estimate in Eq. (15.4.3.1); the radius relation retains its two-phase \(1-\widehat\alpha_e\) denominator and is an approximation in a three-phase mixture.

15.4.5. Low-Pressure Source Enhancement

An optional multiplier increases the source magnitude below a prescribed pressure \(p_{\min}>0\):

(15.4.5.1)\[\begin{split}f_p(p)= \begin{cases} 1,&p>p_{\min}\\[2pt] \displaystyle f_0+(1-f_0) \frac{\exp(\eta p/p_{\min})-1}{\exp(\eta)-1},&p\leq p_{\min} \end{cases}\end{split}\]

Here \(f_0\geq1\) is the factor at zero pressure and \(\eta>0\) controls the shape. Replace \(\dot m_{lv}\) by \(f_p\dot m_{lv}\) consistently in phase, energy, and pressure source terms. The neutral choice is \(f_0=1\).

When low pressure drives evaporation, this modification strengthens the expansion source, but does not guarantee positive pressure.

15.5. Pressure Coupling and Bounded Phase Transport

15.5.1. A Weighted Sum of the Phase Balances

The cavitation pressure equation enforces a weighted combination of the three phase mass balances. The weights express the volume change associated with transferring mass between phases of different density.

For cell \(P\) of volume \(\Omega_P\), let \(M_{i,P}=\rho_PY_{i,P}\Omega_P\) be phase mass. On a fixed mesh, write a backward time derivative as \(\mathcal D_tM_i=c_0M_i^{n+1}-H_i\), with \(H_i=c_1M_i^n-c_2M_i^{n-1}\). Backward Euler has \((c_0,c_1,c_2)=(1,1,0)/\Delta t\); constant-step BDF2 has \((c_0,c_1,c_2)=(3/2,2,1/2)/\Delta t\). Variable-step coefficients are given in Higher-Order Temporal Integration and Splitting.

With a common face mass flux \(\dot m_f\), each phase balance is

(15.5.1.1)\[\mathcal D_tM_{i,P} +\sum_{f\in P}\sigma_{Pf}\dot m_fY_{i,f}=\Omega_PS_{i,P}\]

The stored face orientation is left to right, and \(\sigma_{Pf}\) is the cell incidence sign from Loci Face Orientation and Cell Incidence. The sources are \((S_l,S_v,S_g)=(-\dot m_{lv},\dot m_{lv},0)\).

Divide each phase equation by a positive cell reference density \(\bar\rho_{i,P}\) and sum. These densities are held fixed during a pressure solve. Normally liquid and vapor use saturation densities and gas uses \(p/(R_gT)\) at the current cell state; prescribed constant reference densities are another normalization. Define

(15.5.1.2)\[W_P=\sum_i\frac{Y_{i,P}}{\bar\rho_{i,P}} \qquad W_{Pf}=\sum_i\frac{Y_{i,f}}{\bar\rho_{i,P}} \qquad \Delta v_P=\frac{1}{\bar\rho_{v,P}}-\frac{1}{\bar\rho_{l,P}}\]

The weighted residual at the current nonlinear state is

(15.5.1.3)\[\mathcal R_P= c_0\rho_P\Omega_PW_P-\sum_i\frac{H_{i,P}}{\bar\rho_{i,P}} +\sum_{f\in P}\sigma_{Pf}W_{Pf}\dot m_f^* -\Omega_P\Delta v_P\dot m_{lv,P}\]

Here \(\dot m_f^*\) is the provisional pressure-equation flux. Composition is lagged from the shared phase-transport reconstruction. Each history mass \(M_i^n\) retains its own time-level composition. Using the current \(Y_i\) in every history term would change the discrete balance.

The source contribution \(\Delta v\dot m_{lv}\) is positive for evaporation when \(\bar\rho_l>\bar\rho_v\). It expresses expansion per time. NCG has no source contribution, but remains in both \(W\) and the time history. In general \(W_P\ne1/\rho_P\), because the reference densities need not be the physical EOS densities. Thus this construction weights the individual phase balances; it is not simply continuity divided by mixture density.

15.5.2. Pressure Correction and Source Derivative

Let \(p'\) be the pressure correction. Holding composition and other coefficients fixed, an optional acoustic linearization is \(\rho'\simeq\chi p'/\widetilde a^2\), where \(\chi\) is one when these compressible Jacobian terms are retained and zero when they are omitted. \(\widetilde a\) is the acoustic speed used in the numerical linearization, possibly restricted by numerical bounds. Omitting these terms changes the approximate Jacobian, not the compressible EOS.

For a face oriented from left cell \(L\) to right cell \(R\), the common pressure-based flux correction has the form

(15.5.2.1)\[\dot m_f'=d_f(p_L'-p_R') +\chi Q_f^*\frac{p_{\mathrm{up}}'}{\widetilde a_{\mathrm{up}}^2}\]

Here \(d_f\) is the mass-flux coefficient from momentum interpolation, \(Q_f^*\) is provisional volume flux, and \(\mathrm{up}\) identifies its upwind cell. The momentum and nonorthogonal corrections are developed in Pressure-Correction Equation on Unstructured Grids; they are held known in this linearization.

Write the source change as \(\dot m_{lv}'\simeq S_p p'\) with temperature, composition, densities, geometry, and source multipliers frozen. Linearizing Eq. (15.5.1.3) gives

(15.5.2.2)\[\left(\frac{\chi c_0\Omega_PW_P}{\widetilde a_P^2} -\Omega_P\Delta v_PS_{p,P}\right)p_P' +\sum_{f\in P}\sigma_{Pf}W_{Pf}\dot m_f' =-\mathcal R_P\]

For nonnegative phase availability, increasing pressure reduces evaporation or increases condensation, giving a nonpositive source derivative. Moving it to the left gives the positive diagonal contribution \(-\Omega_P\Delta v_PS_{p,P}\) when \(\Delta v>0\). For Merkle, away from the branch junction,

(15.5.2.3)\[\begin{split}S_p= \begin{cases} -\rho Y_l/(\tau_+t_rq_r),&p<p_s\\ -\rho Y_v/(\tau_-t_rq_r),&p>p_s \end{cases}\end{split}\]

The square-root models use a centered pressure difference,

(15.5.2.4)\[S_p\simeq\frac{\dot m_{lv}(p+\delta p)-\dot m_{lv}(p-\delta p)} {2\delta p}\]

The finite perturbation \(\delta p>0\) regularizes the singular derivative at zero pressure departure. The frozen quantities include NCG density in the fixed-count and Zwart availability factors, and bubble geometry and gas pressure in the transported-count model. If the low-pressure multiplier is used, it multiplies \(S_p\) but is also held fixed. Thus \(S_p\) approximates the source derivative at fixed state properties.

The weighting \(W_{Pf}\) belongs to the pressure-equation row. The actual face mass-flux correction remains the shared Eq. (15.5.2.1); it is not multiplied separately by the two neighboring cell weights. Those weights can differ, making the pressure matrix nonsymmetric without breaking pairwise mass-flux conservation.

15.5.3. Phase and Temperature Updates

The phase predictor treats a destruction term implicitly. With other quantities lagged, write the source as \(\dot m_{lv}=S_C-K_dY_v\) and move \(\Omega_PK_dY_v\) to the diagonal. For Merkle, \(S_C=\rho k_+(1-Y_g)\) and \(K_d=\rho(k_++k_-)\), so this coefficient is nonnegative. Other source models use an analogous split while retaining their lagged geometry and availability factors. The split is a numerical linearization; it does not define a different converged source law.

Within PIMPLE, the usual pointwise phase and temperature correctors reevaluate the mass-transfer rate at corrected pressure and use it explicitly with respect to the new phase fraction and temperature. The same rate enters vapor mass with sign \(+\). In the temperature equation divided by \(c_p\), its factor is \(-L/c_p\). Both backward-Euler and BDF2 pressure-work contributions retain \(\beta_{T,m}T\). The general corrector sequence is given in PIMPLE Predictor–Corrector Method.

An optional local coupled solve is available for Merkle. With neighbors, pressure, density, gas fraction, and transport coefficients held fixed, its two residuals can be written

(15.5.3.1)\[\begin{split}\begin{aligned} r_T(T,Y_v)&={}^{(T)}\!a_PT-{}^{(T)}\!b_P^{\mathrm{loc}} +\frac{\Omega_P}{c_p}L(T)\dot m_{lv}(T,Y_v)\\ r_Y(T,Y_v)&={}^{(Y_v)}\!a_PY_v-{}^{(Y_v)}\!b_P^{\mathrm{loc}} -\Omega_P\dot m_{lv}(T,Y_v) \end{aligned}\end{split}\]

Here \({}^{(T)}\!a_P\) and \({}^{(Y_v)}\!a_P\) are the assembled diagonals excluding phase change. The local right-hand sides \({}^{(T)}\!b_P^{\mathrm{loc}}\) and \({}^{(Y_v)}\!b_P^{\mathrm{loc}}\) include lagged neighbor contributions and exclude phase change. A local Newton step solves the \(2\times2\) Jacobian system for \(\delta T,\delta Y_v\). Temperature changes the tabulated saturation pressure and latent heat; at fixed temperature, \(\partial\dot m_{lv}/\partial Y_v=-\rho(k_++k_-)\) for Merkle. The local solve uses the saturation-table threshold. Pressure, EOS density, and neighboring states are updated by the surrounding flow iteration.

15.5.4. Shared Face Composition

For an admissible upwind cell, let \(\delta Y_v\) and \(\delta Y_g\) be the limited spatial reconstruction increments and set \(\delta Y_l=-\delta Y_v-\delta Y_g\). A common scale preserves all three phase bounds:

(15.5.4.1)\[\theta_f^{\mathrm{rec}}= \min\left(1,\min_{i:\delta Y_i<0} \frac{Y_{i,\mathrm{up}}}{-\delta Y_i}\right) \qquad Y_{i,f}=Y_{i,\mathrm{up}}+\theta_f^{\mathrm{rec}}\delta Y_i\]

An empty inner minimum imposes no restriction. The fractions remain nonnegative and sum to one. First-order transport uses zero increments. Incoming boundary flow uses the prescribed boundary composition; outgoing flow uses the interior reconstruction. The same composition is shared by vapor transport, NCG transport, and pressure weighting at each stage.

These are face bounds. They alone do not ensure that the mass removed through several faces and the source term leave an admissible cell state.

15.5.5. Conservative Flux Rescaling

An optional second restriction limits the total outgoing correction by the available phase inventories. Let \(Y_{i,f}^{H}\) be the bounded high-order face value above and \(Y_{i,f}^{L}\) its donor-cell value. Define the deferred mass-flux correction

(15.5.5.1)\[C_{i,f}=\dot m_f(Y_{i,f}^{H}-Y_{i,f}^{L}) \qquad C_{l,f}=-C_{v,f}-C_{g,f}\]

For each cell, sum only corrections that remove phase mass:

(15.5.5.2)\[L_{i,P}=\sum_{f\in P}[\sigma_{Pf}C_{i,f}]_+\]

Incoming corrections are not counted as available mass because the neighboring cell may limit them. The allowable loss must therefore come from the low-order transport and source update.

For the ordinary PIMPLE corrector on a fixed mesh, define the cell mass \(M_P^n=\rho_P^n\Omega_P\) and total incoming mass flux \(I_P=\sum_{f\in P}[-\sigma_{Pf}\dot m_f]_+\). Subtracting \(Y_i\) times mixture continuity from the phase balance leaves the temporal term \((c_1M_P^n-c_2M_P^{n-1})Y_{i,P}-H_{i,P}\). The pointwise corrector lags the \(c_2\) term, neighbor fractions, and phase source. Its common phase diagonal and low-order numerator are

(15.5.5.3)\[\begin{split}\begin{aligned} {}^{(Y)}\!a_P^{\mathrm{corr}}&=c_1M_P^n+I_P\\ B_{i,P}^{L}&=H_{i,P}+c_2M_P^{n-1}Y_{i,P}^{\mathrm{lag}}\\ &\quad+\sum_{f\in P}[-\sigma_{Pf}\dot m_f]_+ Y_{i,\mathrm{up},f}^{\mathrm{lag}}+\Omega_PS_{i,P}^{\mathrm{lag}}\\ Y_{i,P}^{L}&=\frac{B_{i,P}^{L}}{{}^{(Y)}\!a_P^{\mathrm{corr}}} \end{aligned}\end{split}\]

Here \(\mathrm{lag}\) denotes known values at this corrector stage; the source uses corrected pressure. Boundary inflow supplies its prescribed composition. The \(c_2\) term is absent for backward Euler.

Since the lagged fractions sum to one and \(\sum_iS_i=0\), \(\sum_iB_{i,P}^{L}={}^{(Y)}\!a_P^{\mathrm{corr}}>0\). Thus nonnegative numerators give nonnegative low-order fractions that sum to one. The implicit predictor also requires nonnegative phase history \(H_{i,P}\). A common allowance for both stages is

(15.5.5.4)\[B_{i,P}=\min\left(H_{i,P},[B_{i,P}^{L}]_+\right) \qquad \theta_P=\min\left(1,\min_{i:L_{i,P}>0}\frac{B_{i,P}}{L_{i,P}}\right)\]

For an interior face, take \(\theta_f=\min(\theta_L,\theta_R)\). The final flux is

(15.5.5.5)\[F_{i,f}=\dot m_fY_{i,f}^{L}+\theta_f C_{i,f}\]

Because \(\theta_f\leq\theta_P\), the total correction removing phase mass satisfies

(15.5.5.6)\[\sum_{f\in P}[\sigma_{Pf}\theta_f C_{i,f}]_+ \leq\theta_PL_{i,P}\leq B_{i,P}\leq B_{i,P}^{L} \qquad (B_{i,P}^{L}\geq0)\]

The corrected fraction therefore obeys

(15.5.5.7)\[{}^{(Y)}\!a_P^{\mathrm{corr}}Y_{i,P}^{\mathrm{corr}} =B_{i,P}^{L}-\sum_{f\in P}\sigma_{Pf}\theta_f C_{i,f}\geq0\]

The corrections sum to zero over phases, so the corrected fractions still sum to one. Each is therefore at most one. One shared \(\theta_f\) also preserves \(\sum_iF_{i,f}=\dot m_f\). The same face flux enters neighboring cells with opposite incidence signs, preserving pairwise phase conservation.

Boundary outflow uses the interior cell scale; prescribed inflow has no reconstructed correction. The predictor and corrector recompute allowances for their respective fluxes and pressure-dependent sources.

Conditions for phase bounds

The derivation assumes a fixed mesh, second-order upwind transport, backward Euler or BDF2, and phase equations with the continuity term subtracted and no artificial phase diffusion. Phase bounds require nonnegative history, an admissible low-order update, nonnegative production/destruction coefficients, and a converged implicit predictor. The explicit bound above applies to the ordinary pointwise corrector.

BDF2 does not guarantee \(H_i\geq0\), even when both previous phase masses are nonnegative. A flux restriction also cannot repair a source step that already removes more phase mass than is available. These conditions limit the usable time step even when face fluxes are conservatively limited.

Clipping a cell fraction changes its phase inventory without a compensating neighbor flux. It is therefore not a conservative bounding operation. At fixed density, changing \(Y_i\) by \(\Delta Y_i\) changes phase mass by \(\rho\Omega_P\Delta Y_i\); subsequent EOS recovery can also change mixture density.

15.6. Transported Nuclei and Bubble Growth

15.6.1. Number Balance and Bubble Geometry

An alternative to prescribing \(n_0\) is to transport \(N_b\), the number of bubbles per mixture volume, in \(\mathrm{m^{-3}}\). Its balance is

(15.6.1.1)\[\frac{\partial N_b}{\partial t}+\nabla\cdot(N_b\vec u) =\nabla\cdot(D_N\nabla N_b)+J+bN_b-\frac12KN_b^2\]

Here \(D_N\geq0\) is an optional number diffusivity, \(J\) a volumetric birth rate, \(b\) a binary-breakup frequency in \(\mathrm{s^{-1}}\), and \(K\) a coalescence kernel in \(\mathrm{m^3\,s^{-1}}\). Each binary breakup adds one bubble. Each coalescing pair removes one; the factor \(1/2\) prevents counting the pair twice. These terms change bubble count, not phase mass. The time term transports \(N_b\Omega_P\), and the advective flux is \(N_{b,f}Q_f\) with volume flux \(Q_f\). There is no extra density multiplier.

Use the physical EOS volume fractions for this closure. With \(\alpha_b=\alpha_v+\alpha_g\), the single-size geometry is

(15.6.1.2)\[R_b=\left(\frac{3\alpha_b}{4\pi N_b}\right)^{1/3} \qquad R_e=\max(R_b,R_0)\]

The radius is zero when count or void volume vanishes; \(R_e\) supplies a finite seed scale for growth. The two-phase geometric conversion between count per mixture volume and count per liquid volume is \(N_b=n_0\alpha_l\).

For transported count, define the inception seed by

(15.6.1.3)\[z_b=\frac{4\pi}{3}N_bR_0^3 \qquad Y_0=\min\left(\frac{\rho_v}{\rho}\frac{z_b}{1+z_b},Y_l\right)\]

With a signed radial speed \(u_R\), the availability-weighted mass source is

(15.6.1.4)\[\dot m_{lv}=\frac{3\rho}{R_e}Y_l \left(\max(Y_v,Y_0)[u_R]_+-Y_v[-u_R]_+\right)\]

It is zero for \(N_b=0\). All densities in Eqs. (15.6.1.2)–(15.6.1.4) come from the physical EOS. The factors \(Y_l\) and \(Y_v\) enforce phase availability within this single-size mass-transfer closure.

15.6.2. Reduced Growth with Surface Tension and Gas Pressure

The basic transported-count growth law uses \(u_R=\operatorname{sgn}(p_s-p)\sqrt{2|p_s-p|/(3\rho_l)}\) with the current liquid EOS density. The reduced Rayleigh–Plesset extension retains viscosity, surface tension, and gas pressure while omitting radial acceleration. The gas pressure supplied by the resolved NCG inventory is

(15.6.2.1)\[p_b=\frac{\rho Y_gR_gT}{N_b(4\pi R_e^3/3)}\]

It is zero when \(N_b=0\). This assumes gas at the mixture temperature occupies the modeled bubbly volume.

Define two dimensional coefficients,

(15.6.2.2)\[A_R=\frac{4\mu_l}{3\rho_lR_e} \qquad Q_R=\frac{2(p_s+p_b-p-2\sigma/R_e)}{3\rho_l} \qquad \Delta_R=A_R^2+Q_R\]

Deleting \(R_b\ddot R_b\) from Eq. (15.4.2.1) gives \(u_R^2+2A_Ru_R-Q_R=0\). The root continuous through zero growth is extended to collapse with the signed-root prescription of Pardo Vigil and coauthors [PVSuarezFernandezGonzalezPerezP25]:

(15.6.2.3)\[\begin{split}u_R= \begin{cases} \displaystyle\frac{Q_R}{\sqrt{\Delta_R}+A_R},&\Delta_R\geq0\\[6pt] -A_R-\sqrt{-\Delta_R},&\Delta_R<0 \end{cases}\end{split}\]

The first expression equals \(-A_R+\sqrt{\Delta_R}\), with less cancellation near zero growth. The second is a continuation of the model, not a real root of the quadratic. This closure changes the growth rate in Eq. (15.6.1.4) while retaining the bulk mixture EOS. It omits radial acceleration and a separate bubble gas-energy history, so it cannot represent collapse rebound.

15.6.3. Population Sources

The birth rate \(J\) can combine homogeneous nucleation within liquid and nucleation at adjacent walls. Both depend on the energy needed to form a vapor nucleus; the wall rate also depends on contact angle. Breakup and coalescence can use prescribed nonnegative \(b,K\) or turbulence-dependent rates based on dissipation and bubble size. The detailed closures are given in Cavitation Population Source Models.

These sources change the mean bubble size through Eq. (15.6.1.2). The single-size approximation does not retain separate daughter sizes or a bubble-size distribution.

15.7. Turbulence Corrections

15.7.1. Density-Function Correction

The density-function correction (DFR) reduces the SST RANS eddy viscosity in a liquid–vapor mixture. Let \(\mu_t^{(0)}\) be the uncorrected eddy viscosity. For two phases, the Reboud-type effective density replaces the linear liquid-volume dependence by a power, as used in cavitating RANS models [GE20]:

(15.7.1.1)\[\rho_D=\rho_v+\alpha_l^{n_D}(\rho_l-\rho_v) \qquad n_D\geq1\]

The NCG extension retains the gas density contribution and applies the power only within the condensable volume. Define \(\alpha_c=\alpha_l+\alpha_v\) and \(q=\alpha_l/\alpha_c\). Then

(15.7.1.2)\[\rho_D=\alpha_c\left[\rho_v+q^{n_D}(\rho_l-\rho_v)\right] +\alpha_g\rho_g\]

Subtracting the physical density \(\rho=\alpha_c[\rho_v+q(\rho_l-\rho_v)]+\alpha_g\rho_g\) and using \(\alpha_l\rho_l/\rho=Y_l\) gives

(15.7.1.3)\[F_\rho=\frac{\rho_D}{\rho} =1-Y_l\left(1-\frac{\rho_v}{\rho_l}\right) \left(1-q^{n_D-1}\right) \qquad \mu_t=F_\rho\mu_t^{(0)}\]

All densities are physical EOS densities at the current \(p,T\). For \(\rho_l>\rho_v>0\) and admissible composition, \(0<F_\rho\leq1\). The factor is one for \(n_D=1\), pure phases, or a mixture with no vapor or no liquid. The uncorrected value is also used outside the subcritical range or when \(\rho_l\leq\rho_v\).

The conventional exponent is \(n_D=10\). It is empirical, and the NCG extension is a modeling choice. DFR changes turbulent transport and its associated production terms, while the EOS and conserved mixture density remain those in Eq. (15.2.1.2). Its scope here is the SST 2003 RANS closure described in RANS-Based Turbulence Models.

15.7.2. Turbulent Saturation-Pressure Correction

A separate optional correction represents unresolved pressure fluctuations by increasing the source threshold:

(15.7.2.1)\[p_s=p_{\mathrm{sat}}(T)+0.195\rho k\]

The quantity \(\rho k\) has pressure units. This correction changes the pressure departure used by the mass-transfer law; it does not change the thermodynamic saturation table. It is independent of the eddy-viscosity correction above.

15.8. Applicability and Limitations

The model describes cavitation at the mixture scale. Its shared velocity and temperature and single-size bubble closures leave phase slip, bubble-scale thermal layers, collapse rebound, and erosion unresolved. Predictions depend on the chosen rate coefficients, nuclei model, turbulence corrections, and valid thermodynamic property coverage.

An accepted time level must satisfy the phase, pressure, thermal, and EOS relations together. Convergence of the pressure correction alone does not establish that consistency; source stiffness and phase-history admissibility can impose additional time-step restrictions.

15.9. References

[Bre95]

Christopher E. Brennen. Cavitation and Bubble Dynamics. Oxford University Press, 1995. URL: https://media.library.caltech.edu/CaltechBOOK:1995.001/.

[GE20]

Linlin Geng and Xavier Escaler. Assessment of RANS turbulence models and Zwart cavitation model empirical coefficients for the simulation of unsteady cloud cavitation. Engineering Applications of Computational Fluid Mechanics, 14(1):151–167, 2020. doi:10.1080/19942060.2019.1694996.

[MFB98]

C. L. Merkle, J. Feng, and P. E. O. Buelow. Computational modeling of the dynamics of sheet cavitation. In Proceedings of the 3rd International Symposium on Cavitation. Grenoble, France, 1998.

[PVSuarezFernandezGonzalezPerezP25]

Álvaro Pardo Vigil, Laura Suárez Fernández, José González Pérez, and Adrián Pandal. An improved, Rayleigh–Plesset based homogeneous cavitation model accounting for microbubble behaviour and turbulent interaction. International Journal of Multiphase Flow, 191:105334, 2025. doi:10.1016/j.ijmultiphaseflow.2025.105334.

[SS01]

Günter H. Schnerr and Jürgen Sauer. Physical and numerical modeling of unsteady cavitation dynamics. In Proceedings of the 4th International Conference on Multiphase Flow. New Orleans, USA, 2001. URL: https://publikationen.bibliothek.kit.edu/20722001.

[ZGB04]

Philip J. Zwart, Andrew G. Gerber, and Thabet Belamri. A two-phase flow model for predicting cavitation dynamics. In Proceedings of the 5th International Conference on Multiphase Flow, number 152. Yokohama, Japan, 2004. URL: https://www.researchgate.net/publication/306205415_A_two-phase_flow_model_for_predicting_cavitation_dynamics.