Simulation of Cavitating Flows

Load the cavitation module in the run control file:

loadModule: cavitation_nist

The module uses tabulated REFPROP data for fluid properties. Liquid, vapor, and non-condensable gas (NCG) share pressure, temperature, and velocity. Vapor and NCG mass fractions are transported separately; this chapter uses \(\alpha_V\) and \(\alpha_{NCG}\) for these mass fractions.

Cavitation Models

Select Merkle, Zwart, or SauerSchnerr with source in cavitationEquationOptions. Singhal is not implemented and is rejected during setup. Enter model parameters as numbers in the units listed below. Example values illustrate syntax; they are not defaults or a calibration. See Cavitating Flows in the Theory Guide for the source equations.

Merkle Source Term Model

The Merkle model [MeFB1998] [HoAU2007] makes phase-change rates proportional to the pressure departure from the local vapor pressure. All five parameters in MerkleSourceParameters are required:

Merkle parameters

Option	Meaning	Units
tauVap	Vaporization time-scale parameter	Dimensionless
tauCond	Condensation time-scale parameter	Dimensionless
L	Reference length	m
V	Reference speed	m/s
rho	Reference density	kg/m3

Larger tauVap and tauCond weaken their respective rates. They are reciprocals of multiplicative rate coefficients: \(\tau=1/C\). For a pressure departure comparable to the reference dynamic pressure, the phase-change time is approximately \(\tau L/V\). Compare parameters across cases using the reference scales as well as \(\tau\).

cavitationEquationOptions: <source=Merkle>
MerkleSourceParameters: <tauVap=1.0, tauCond=80.0, L=0.7, V=1.5, rho=1.15>

Zwart Source Term Model

The Zwart model [GeZB2004] uses a prescribed bubble radius and nuclei volume fraction, with separate vaporization and condensation rate multipliers. All four parameters in ZwartSourceParameters are required:

Zwart parameters

Option	Meaning	Units
fVap	Vaporization rate multiplier	Dimensionless
fCond	Condensation rate multiplier	Dimensionless
RB	Bubble radius	m
aNuc	Nuclei volume fraction	Dimensionless

cavitationEquationOptions: <source=Zwart>
ZwartSourceParameters: <RB=1.0e-6, aNuc=5.0e-4, fVap=50.0, fCond=0.01>

Sauer-Schnerr Source Term Model

The Sauer-Schnerr model [ScSS2008] derives the bubble radius from vapor content and nuclei density. For fixed nuclei density, use SauerSchnerrSourceParameters:

Sauer-Schnerr parameters

Option	Meaning	Units	Default
rNuc	Seed radius used to initiate vaporization; finite and positive	m	1.0e-5
n	Prescribed nuclei number density; finite and positive	m-3	Required for fixed nuclei density

cavitationEquationOptions: <source=SauerSchnerr>
SauerSchnerrSourceParameters: <rNuc=1.0e-5, n=1.0e8>

Transported nuclei density uses nb instead of fixed n; see Transported Nuclei and Bubble Population.

Zwart and fixed-nuclei Sauer-Schnerr use square-root pressure dependence. The signed bubble speed is positive for vaporization and negative for condensation:

\[\begin{split}u_b = \begin{cases} \sqrt{\dfrac{2(P_{vap}-P)}{3\rho_L}}, & P<P_{vap} \\ 0, & P=P_{vap} \\ -\sqrt{\dfrac{2(P-P_{vap})}{3\rho_L}}, & P>P_{vap} \end{cases}\end{split}\]

Here \(P\) is local pressure, \(P_{vap}\) is vapor pressure at the local temperature, and \(\rho_L\) is saturation liquid density. The corresponding vapor source has the same sign.

Source Term Scaling Factor

Below cavitationMinPressure, the source is multiplied by

\[\begin{split}F(P) = \begin{cases} F_{\max} + (1-F_{\max})\dfrac{\exp(kP/P_{\min})-1}{\exp(k)-1}, & P\leq P_{\min} \\ 1, & P>P_{\min} \end{cases}\end{split}\]

Source scaling parameters

Option	Meaning	Default
cavitationMaxSourceFactor	Dimensionless factor \(F_{\max}\) at zero pressure	1
cavitationSourceFactorExponent	Dimensionless exponent \(k\); must be nonzero	1
cavitationMinPressure	Activation pressure \(P_{\min}\) in Pa; must be nonzero	1

The default source factor of one disables scaling. Values greater than one amplify the source below the activation pressure; values between zero and one damp it. cavitationMinPressure sets the activation threshold, not a pressure floor.

Control File Setup

Cavitation requires compressible mode and form=temperature in energyEquationOptions. Total enthalpy and total energy are not supported. Provide vapor_y and NCG_y at dedicated inlets and in the initial condition. Use NCG_y=0 when no non-condensable gas is present.

Boundary Conditions

Dedicated inlets are subsonicInlet, supersonicInlet, and inflow. For example:

Inlet=subsonicInlet(v=10.11 m/s, T=300.0 K, k=0.0, omega=1000.0, vapor_y=0.0, NCG_y=1.0e-06)

For outlet boundaries, the normal Mach number should be evaluated using the local mixture sound speed, which can decrease substantially as vapor forms. The fixedPressureOutlet condition should be used only where the normal flow remains subsonic.

The outflow and farfield conditions recover boundary temperature from the complete liquid/vapor/NCG mixture EOS:

\[\frac{1}{\rho} = \frac{1-\alpha_V-Y_g}{\rho_L(p,T)} + \frac{\alpha_V}{\rho_V(p,T)} + Y_g\frac{R_g T}{p}.\]

Here \(Y_g\) is the NCG mass fraction and \(R_g\) follows the configured NCG molecular weight. The inversion uses the pressure assigned to the face and searches within TclipMin and TclipMax. If a bracket cannot be found or the EOS residual does not converge, Stream stops with a boundary inversion diagnostic. It does not accept an unconverged temperature and then clip it. The property tables retain their existing endpoint-clamping behavior; these temperature bounds do not extend the tables’ physical validity.

The two boundaries handle incoming flow differently:

• outflow retains its existing blocked-flow branch: when the adjacent cell’s normal velocity points inward, the face velocity is set to zero, composition is extrapolated, and the face pressure is the smaller of the cell pressure and the specified pressure. This is not a prescribed incoming reservoir state. Use farfield when an incoming reservoir must be defined.

• For subsonic farfield, the sign of the updated face-normal velocity selects the advected state. Outgoing fluid uses the interior composition and thermal state. Incoming fluid uses the specified reference composition, reference thermal state, and tangential velocity. Its density is corrected acoustically from the reference pressure to the face pressure before EOS inversion, so its face temperature need not equal the specified reference temperature. Supersonic outflow extrapolates the interior state; supersonic inflow assigns the reference state.

For the velocity form of farfield, specify p, T, v, and vapor_y. Its optional NCG_y defaults to zero; specifying it explicitly is recommended for clarity. For example:

Open=farfield(p=100000 Pa, T=300 K, v=[1 m/s,0,0], vapor_y=0, NCG_y=1e-6)

fixedPressureOutlet does not use this inversion. It continues to extrapolate composition during reversed flow and does not prescribe a separate incoming reservoir composition.

Initial Conditions

Specify the initial composition with the flow state:

initialCondition: <vapor_y=0, NCG_y=0.0, v=10.11 m/s, p=28515.8 Pa, T=300.0 K, k=0.0, omega=1000>

For ql, qr, and initialConditionRegions, see Appendix: Initial Conditions.

Equation of State and Transport Properties

Use cav_tab_tool to generate REFPROP liquid (.liq), vapor (.vap), and saturation (.sat) tables. Set liquid_model and vapor_model to their filename prefixes; the saturation file uses the liquid prefix. transport_model: module selects transport properties from these tables.

liquid_model: NITROGEN
vapor_model: NITROGEN
transport_model: module

By default, the tabulator rejects ranges that cross the phase boundary. --extend_to_phase_boundary supplies a numerical continuation for rectangular tables; it does not create physical pure-phase states beyond that boundary. See Cavitation REFPROP Tabulation Utility (cav_tab_tool) for usage and limits.

Cavitation Equation Options

Add solver options to the selected source model:

cavitationEquationOptions: <linearSolver=SGS, relaxationFactor=0.9, maxIterations=5,
                            source=Zwart, vaporProductionRateRelaxationFactor=0.7>

Some Cavitation Module Options

Parameter	Description
AlphaViscosity	Artificial viscosity added to vapor and NCG mass fraction transport equations for numerical stability purposes.
cavitationInviscidFlux	Flux scheme used for vapor and NCG transport equations. Supported values are: FOU and SOU.
liquid_rho	Reference value for liquid-phase density used in normalizing the liquid term in the pressure-correction equation.
vapor_rho	Reference value for vapor-phase density used in normalizing the vapor term in the pressure-correction equation.
NCG_rho	Reference value for non-condensable-gas density used in normalizing the NCG term in the pressure-correction equation.
turbulentVapourPressureCorrection	Activates the turbulent correction to the vapor pressure.

See Appendix: Linear Solvers for solver settings; the default is SGS. Set cavitation: off to disable vapor source terms (default on). Omitted reference densities, or values of -1, use local phase densities to normalize the pressure-correction equation.

Phase Transport and Convergence

Vapor and NCG convection use the same face composition as the cavitation pressure-correction equation. With cavitationInviscidFlux: SOU, Stream reconstructs both mass fractions from the upwind cell and scales their increments by a common factor when needed to keep

\[\alpha_V \geq 0,\qquad \alpha_{NCG} \geq 0,\qquad \alpha_V+\alpha_{NCG}\leq 1.\]

The donor is selected using volume flux relative to the moving grid. Outgoing boundary flow uses the interior reconstruction; incoming flow uses the composition supplied by that boundary condition. This does not add a reservoir state to boundaries that extrapolate composition during backflow. The boundary EOS and the pressure equation’s reference-density normalization are unchanged. PIMPLE pressure correction uses lagged composition; its phase corrector selects the donor again using the corrected relative flux.

This face limiter preserves admissible face fractions. It does not guarantee an admissible cell update. The existing cell limits can still prevent a sharp second-order solution from satisfying the discrete phase equations. Such a case can now fail the convergence check described below.

Set conservativeFluxRescale=on in cavitationEquationOptions to also limit the conservative higher-order flux correction using cell inventories. The default is off. One shared face factor protects vapor, NCG, and the remaining liquid; internal corrections enter neighboring cells with opposite signs. Inflow retains the boundary composition, and outflow corrections use the interior cell’s allowance.

cavitationEquationOptions: <source=Merkle,linearSolver=SGS,
    relaxationFactor=0.3,maxIterations=20,conservativeFluxRescale=on>
cavitationInviscidFlux: SOU
timeIntegrator: BDF
PDEContinuityAdded: yes
AlphaViscosity: 0.0

This option supports fixed meshes, BDF or BDF2, and the Merkle, Zwart, and Sauer–Schnerr sources, with phase change enabled or disabled. Sauer–Schnerr may use a prescribed nuclei density or transported nb, including its bubble-growth options. It supports SIMPLE/SIMPLEC and the ordinary PIMPLE phase corrector. SOU, PDEContinuityAdded=yes, and AlphaViscosity=0 are required. Moving meshes, other time integrators, and explicitCorrectorNewtonIteration=on are rejected. The reconstruction limiter selected by limiter still applies.

BDF2 additionally requires a nonnegative weighted history for each constituent. For equal timesteps, this requires 4 M_old - M_older >= 0 for each cell’s vapor, NCG, and liquid mass. A material violation stops the run before the history is used. Reduce the timestep or restart with BDF. Rescaling does not clip that history or switch the time integrator automatically.

Rescaling reduces the higher-order correction where an inventory is small; it can increase numerical diffusion. It does not change the EOS or cavitation source model. PIMPLE’s explicit source update must still stay within the available phase inventory. A material bound violation in that corrector stops the run with a diagnostic rather than clipping the result. Reduce the timestep if the explicit source step is too large. Phase-equation convergence checks remain active with either value of conservativeFluxRescale.

Add cavitationFlux to diagnostics to inspect the predictor:

diagnostics: <cavitationFlux>
plot_output: vapor_y,NCG_y,phaseFluxReduction

The CFLUX: line reports cycle, iteration, cell count, limited-cell count, mean correction reduction, material-clipped cell count, and absolute vapor and NCG clipping adjustments in kg. The material-clipped count excludes changes no larger than 64 machine epsilons in either mass fraction. It follows print_freq. These are per-iteration measurements; summing the clipping adjustments over iterations does not measure physical mass loss. phaseFluxReduction ranges from zero (unrestricted) to one (no higher-order correction). Its plotted nodal value is the largest adjacent cell value. The diagnostic also works with rescaling off, when its correction reduction is zero.

Both phase equations participate in nonlinear convergence, including when cavitation: off disables vapor production. Optional vapor and NCG entries override default in the usual tolerance dictionaries. For example:

convergenceAbsoluteTolerance: <default=1e-10,vapor=1e-12,NCG=1e-14>
convergenceRelativeTolerance: <default=1e-6,vapor=1e-8,NCG=1e-8>

These are illustrative values, not universal tolerances. The absolute phase residual is the sum of the absolute cell equation imbalances, in kg/s. Each phase passes when either its absolute tolerance or its relative tolerance is met. The relative reference is the first positive phase residual within the current physical timestep; it is reset at each timestep, including restart entry. An initially zero residual does not cause division by zero. A small roundoff allowance, based on the assembled equation terms, also permits balances that cannot be reduced further in floating-point arithmetic.

An entry omitted from a dictionary inherits its default value; if neither is present, its value is 1e-30. The legacy convergenceTolerance retains its existing precedence: an active legacy value supplies the absolute tolerance and sets the relative tolerance to 1e-30. Phase tolerances must be finite and nonnegative.

For PIMPLE, the phase stopping decision uses the state after the final inner correction. Stream also checks the final state before accepting the timestep and writing its phase restart data. If either phase remains unconverged at maxIterationsPerTimeStep, the run stops with the phase name, timestep, residual, and tolerance. A final SIMPLE update that invalidates an earlier stopping decision is also rejected. The solver does not automatically retry the failed timestep or reduce its size. Check the residual history, phase bounds, timestep, and spatial discretization when diagnosing either failure.

Transported Nuclei and Bubble Population

The optional nuclei equation replaces the prescribed Sauer–Schnerr bubble count with a transported field. Here nb is the number of nuclei per cubic metre of mixture. It is not a mass fraction. Enable it with:

SauerSchnerrSourceParameters: <rNuc=3e-6, nucleationTransport>
nucleationEquationOptions: <linearSolver=SGS, relaxationFactor=0.5,
                            maxIterations=20, nucleation=none,
                            population=none, bubbleGrowth=rayleigh>
nucleationInviscidFlux: FOU

Specify a finite, nonnegative nb in the initial condition and at each inlet. Do not also specify the constant n. The fixed-n model retains its existing reference-volume convention; its input is not interchangeable with transported nb. BDF and BDF2 are supported. Start with FOU; SOU is available, but a materially negative number density stops the run and requires a smaller timestep or a more bounded discretization.

Writing \(N=\mathtt{nb}\), the number equation is

\[\frac{\partial N}{\partial t}+\nabla\cdot(N\mathbf{u}) =\nabla\cdot(D_N\nabla N)+J+bN-\frac12 K N^2.\]

The finite-volume equation uses relative volume flux on a moving mesh and stores \(NV\), including both histories for BDF2. nbViscosity is \(D_N\) in \(\mathrm{m^2/s}\) and defaults to zero. The nonlinear stopping check includes the number residual. absoluteTolerance in nucleationEquationOptions defaults to \(10^{-3}\) nuclei/s; relativeTolerance defaults to \(10^{-8}\). A roundoff allowance is scaled by the assembled number equation. An unconverged number equation at the nonlinear iteration limit stops the run.

Outlets extrapolate nb unless an incoming value is supplied. An outlet nb then specifies the reservoir state during backflow. Farfield boundaries require a supplied nb. At impermeable walls the advective number flux is zero.

Bubble size and phase change

The local EOS gives the liquid and vapor densities at the current pressure and temperature. NCG uses its ideal-gas EOS. These densities and the transported mass fractions give the bubbly volume fraction \(\phi_B=\phi_V+\phi_G\). The representative radius is

\[R=\left(\frac{3\phi_B}{4\pi N}\right)^{1/3},\qquad R_e=\max(R,r_{nuc}).\]

The rate calculation uses \(R_e\). Zero nuclei gives zero phase-change rate. rNuc supplies a seed size; it does not add vapor or gas mass. The transported model uses nonnegative evaporation and condensation coefficients in the vapor equation and uses the same net mass source in pressure correction and latent heat. Gas is never converted into liquid.

bubbleGrowth=rayleigh is the default. It uses the usual inertial square-root growth/collapse law. bubbleGrowth=reducedRP adds surface tension, liquid viscosity, and isothermal gas pressure to the algebraic Rayleigh–Plesset approximation:

\[ \begin{align}\begin{aligned}a=\frac{4\mu_L}{3\rho_L R_e},\qquad q=\frac{2}{3\rho_L}\left(p_V+p_G-p-\frac{2\sigma}{R_e}\right), \qquad \chi=a^2+q,\\\dot R=-a+\operatorname{sgn}(\chi)\sqrt{|\chi|}.\end{aligned}\end{align} \]

The positive-discriminant branch is evaluated in a form that avoids cancellation near equilibrium. Surface tension comes from the saturation table and viscosity from the liquid EOS table. The gas mass per bubble is \(m_G=\rho Y_G/N\), so \(p_G=m_G R_G T/(4\pi R_e^3/3)\); the density in this expression is reconstructed from the current EOS state.

This is an experimental Stream adaptation of the reduced growth law in Pardo Vigil et al. (2025). It retains Stream’s three-component mixture EOS and common temperature. It does not implement the paper’s separate polytropic bubbly-phase EOS, slip model, or turbulence-pressure correlation. In particular, it cannot reproduce all gas-supported microbubble equilibria of that model. Radial inertia, collapse rebound, and acoustic emission are not resolved. Use measured nuclei information and case-specific validation before predictive use.

Population sources

population=none transports existing nuclei without breakup or coalescence. population=constant accepts breakupRate in \(\mathrm{s^{-1}}\) and coalescenceKernel in \(\mathrm{m^3/s}\), both defaulting to zero. This also provides a simple reference model for a measured or prescribed rate.

population=turbulentSingleSize uses a single-size reduction of the Prince–Blanch turbulent collision model and an equal-volume binary Luo–Svendsen breakup approximation. It currently requires menterSST2003 RANS. Dissipation is \(\epsilon=\beta^*k\omega\). breakupMultiplier and coalescenceMultiplier default to one; filmThicknessRatio defaults to \(10^4\) and must exceed one. Setting either multiplier to zero disables that contribution.

Breakup adds one bubble per binary event. Coalescence removes one bubble per unordered pair, which gives the factor \(1/2\) in the number source. Both change the representative radius without adding phase mass. This model does not retain a size distribution, daughter sizes, slip, or bubble histories. The breakup expression is a reduced closure, not the full Luo–Svendsen population balance. Rates require validation for the liquid and flow of interest. Large growth rates can require smaller timesteps even when the fluid equations converge.

nucleation=classical adds the inherited classical homogeneous and wall nucleation laws. It requires SpeciesMolecularMass in kg/kmol. The default is nucleation=none. Homogeneous rates are weighted by liquid volume; wall rates by wetted area and the adjacent liquid fraction. The bubble contact angle thetaB defaults to 170 degrees and must lie in \((0,180]\). nbSourceClipMax limits the cell-averaged birth rate in nuclei/\(\mathrm{m^3/s}\) and defaults to \(10^{30}\). These birth laws are highly sensitive to pressure and wetting and are not a substitute for measured inlet nuclei.

The repaired transport uses nbMixture_hdf5 restart files and, for BDF2, nbMixturePrev_hdf5. Older nb_hdf5 files are rejected because the older density-weighted equation had a different conservation convention. Reinitialize nuclei when migrating an older calculation.

DFR Turbulence Correction

The optional Reboud density correction reduces turbulent eddy viscosity in liquid-vapor mixtures. It is available in cavitation_nist with flowRegime: turbulent and model=menterSST2003 RANS. Select it with turbulenceCorrection=DFR in cavitationEquationOptions:

turbulenceEquationOptions: <model=menterSST2003, linearSolver=SGS,
                            relaxationFactor=0.5, maxIterations=20>
cavitationEquationOptions: <source=Merkle, turbulenceCorrection=DFR,
                            linearSolver=SGS, relaxationFactor=0.3,
                            maxIterations=20>
DFRTurbulenceCorrectionParameters: <n=10>

Omitting turbulenceCorrection leaves the ordinary turbulence model unchanged. The exponent n defaults to 10 and must be finite and at least 1. Setting n=1 selects the original eddy-viscosity rules. Other turbulence models and the des, DESLimiter, and rhok options are rejected when DFR is selected. DFR can be used with any supported cavitation source model and also with cavitation: off to study a prescribed mixture.

For a mixture without NCG, the effective density is

\[\rho_{DFR}=\rho_V+\phi_L^n(\rho_L-\rho_V), \qquad F_\rho=\frac{\rho_{DFR}}{\rho}.\]

Here \(\phi_L\) is the liquid volume fraction, whereas Alpha in the code is the vapor mass fraction. The phase densities are evaluated at the local pressure and temperature using the same tables as the mixture EOS. The SST eddy-viscosity rules use \(F_\rho\rho\) in place of \(\rho\), before applying the existing viscosity cap. The physical mixture density, EOS, and mass and momentum storage terms are unchanged.

For NCG, Stream applies the correction within the liquid-vapor portion and retains the gas density contribution. This is a Stream extension of the two-phase correction:

\[ \begin{align}\begin{aligned}\phi_c=\phi_L+\phi_V, \qquad q=\frac{\phi_L}{\phi_c},\\\rho_{DFR}=\phi_c\left[\rho_V+q^n(\rho_L-\rho_V)\right] +\phi_{NCG}\rho_{NCG}.\end{aligned}\end{align} \]

The volume fractions satisfy \(\phi_i=\rho\alpha_i/\rho_i\), where \(\alpha_i\) denotes the phase mass fraction used elsewhere in this chapter. The factor is one when either liquid or vapor is absent, including pure NCG. It is also one at or above the upper saturation-table temperature or when the EOS gives \(\rho_L\leq\rho_V\). The density-contrast limit is continuous and needs no division by \(\rho_L-\rho_V\).

Cell and boundary-face factors use their respective local states. Wall-function cell eddy viscosity receives the same factor before its existing cap; the wall shear law and wall values of \(k\) and \(\omega\) retain their existing definitions. Explicitly prescribed boundary eddy viscosity (a profile or muRatio) remains a user override. The existing rhoFactor plot variable shows the cell factor interpolated to nodes.

turbulentVapourPressureCorrection is a separate option that changes the phase-change pressure threshold. DFR does not enable it. Reducing eddy viscosity can allow cavity shedding, so a steady result is not guaranteed. The exponent is empirical; its default is a starting value, not a universal calibration. See the published SST cavitation assessment for an application of the Reboud correction.

Non-Condensable Gas Specification

All cavitation simulations involve NCG terms in the governing equations, even if the user does not use NCG. The default fluid properties used for the nonCondensableGasProperties variable is air as shown below. To use a different fluid, the user must specify its molar mass m, specific heat at constant pressure cp, laminar dynamic viscosity mu, and thermal conductivity kcond in the nonCondensableGasProperties variable.

nonCondensableGasProperties: <m=28.9647 g/mol, cp=1005 J/kg/K, mu=1.805e-5 kg/m/s,
                              kcond=0.02476 W/m/K>

Output Variables

The table below shows additional field variables that are available for output in cavitating flow simulations. Default variables are output automatically and do not need to be specified in the plot_output line in the run control file.

Field Variable Output for Cavitating Flows.

Variable	Description	Default
bt	Bulk expansion coefficient * temperature	No
F	Source term scaling factor	No
liquid_vof	Liquid volume fraction	No
liquid_y	Liquid mass fraction	No
NCG_vof	NCG volume fraction	No
NCG_y	NCG mass fraction	No
NCG_yResidual	NCG mass fraction eq. residual	No
pV	Vapor pressure	No
rhoL	Liquid density	No
rhoV	Vapor density	No
vapor_y	Vapor mass fraction	No
vapor_vof	Vapor volume fraction	No
vapor_yResidual	Vapor mass fraction eq. residual	No

Diagnostic Units and Interpretation

AlphaProductionRate is an optional signed field in kg/s. It is the volume-integrated vapor source evaluated at the current pressure: positive for evaporation and negative for condensation. AlphaProductionTimeScale is an optional nonnegative field in seconds. With \(S\) denoting that source, it reports the local time to consume the available phase at a frozen source rate:

\[\begin{split}\tau = \begin{cases} \rho V(1-\alpha_V-Y_g)/S, & S>0,\\ \rho V\alpha_V/(-S), & S<0. \end{cases}\end{split}\]

For zero source, including cavitation: off, the time is reported as \(10^{30}\) s. Larger times are capped at this value. This diagnostic does not prescribe a stable timestep or predict the nonlinear evolution of the source. Both fields use the existing maximum-magnitude cell-to-node output interpolation; the formula above describes the cell value.

Output change: AlphaProductionTimeScale now contains seconds. Older versions attempted to encode direction in a signed logarithm and omitted cell volume. Do not interpret older files with the new definition; obtain direction from AlphaProductionRate in newly generated output.

For integrated boundary output, mixture specific internal energy is

\[e_f=(1-\alpha_V-Y_g)e_L(p_f,T_f)+\alpha_V e_V(p_f,T_f) +Y_g(c_{p,g}-R_g)T_f.\]

The liquid/vapor tables supply internal energy directly. Their energy datums are retained; the calorically perfect NCG uses \(e_g(0)=0\). The energy column in flux_<boundary>.dat includes convected internal and kinetic energy using velocity relative to the grid, plus the existing heat transfer contribution. It is not a complete total-energy balance: pressure work and viscous work are not added by this diagnostic.

The component mass diagnostic uses the same relative-grid flux convention and differences \(\rho V\) at matching history levels. Its signed defect is \(D_t\sum\rho V+\sum\dot m_{out}\), in kg/s, using the predictor flux volFluxStar. It should approach zero as continuity converges. BDF2 uses the configured variable-step history coefficients. This diagnostic is not a replacement for an independent accepted-state mass budget.

Cavitation REFPROP Tabulation Utility (cav_tab_tool)

The Python REFPROP tabulation utility is located under in the bin/ directory in the source code. The name of the utility is cav_tab_tool. This tool is for building the liquid, vapor, and saturation input files used by the cavitation_nist workflow.

This is not the same workflow as the tabulated_data refprop utility. Use cav_tab_tool when a cavitation case needs .liq, .vap, and .sat files. Use tabulated_data when a tabular-EoS case needs an EoSTables/<table_name>/ directory.

The Python tabulation utility has a functionality for making calls directly to REFPROP 11. This allows for the tool to tabulate as many pressure levels as the user desires. The user specifies the REFPROP material name, the pressure range, temperature range, and the number of data points to use over each range. The user can also pass a flag (--log_p) to uniformly distribute the pressure points in the log-space to ensure that an equal number of points resolve each pressure decade (this is good for large pressure ranges). The tool also supports adaptive temperature spacing (--adaptive_t) that places more points where the data is highly curved and fewer points where the data is nearly linear.

A set of .liq, .vap, and .sat files are generated by the tool, and these are the inputs used by the cavitation_nist module in Stream. Sample calls to the utility for the liquid/vapor files and the saturation files are given below.

There are two modes that the tool can be used in. One is for tabulating properties outside of a two-phase region in the thermodynamic space. The other is for tabulating the saturation curve properties of the liquid and vapor states.

Liquid/Vapor Tabulation

To generate tabulated liquid and vapor files, one should pass the liquidVapor argument to the tabulation tool. The presence of this argument enables a required set of arguments which are described below:

Liquid/Vapor Tabulation Arguments

Argument	Description
t_min	The minimum temperature range value
t_max	The maximum temperature range value
num_t	The number of temperature levels
p_min	The minimum pressure range value
p_max	The maximum pressure range value
num_p	The number of pressure levels
log_p	An optional argument that distributes the pressure levels uniformly across the logarithmic space of the pressure range
adaptive_t	Optional flag to enable adaptive temperature spacing
extend_to_phase_boundary	Optional flag to write a single phase-boundary continuation row when the requested temperature range contains no stable pure-phase state at a pressure level
t_rel_err	Relative error scale for adaptive spacing
t_abs_err	Absolute error scale for adaptive spacing
t_integrated_err	Target normalized integrated error used to stop adaptive refinement
t_min_points	Minimum number of temperature points used by adaptive spacing
t_max_points	Maximum number of temperature points used by adaptive spacing

If the log_p argument is not provided, a linear tabulation in the pressure dimension is performed. The logarithmic spacing enabled by the log_p argument creates a uniform spacing in the logarithmic space. This type of spacing is useful when there is a large variation in the pressure range, as a linear tabulation would result in a low resolution in the low-pressure region of the tabulation. An example command to generate liquid and vapor tabulation files for nitrogen in the temperature range of 65K-300K and pressure range of 50kPa-70MPa is shown below.

<StreamInstallDir>/bin/cav_tab_tool liquidVapor --fluid_name NITROGEN --t_min 65 --t_max 300 --p_min 50000 --p_max 70e6 --log_p --num_p 200 --num_t 200

In the example above, the liquidVapor argument sets the tool into the mode for tabulating .liq and .vap tables. The fluid_name argument must match an available REFPROP FLD file in the REFPROP fluid database. The tabulation levels are done at constant pressure, and so the t_min and t_max arguments are the temperature range over which to tabulate properties. The p_min and p_max arguments are the inclusive bounds of the pressure range over which you want to tabulate data. The log_p argument converts the pressure range into a logarithmic range and distributes points equally in that space and then converts back to pressure levels. This ensures an equal number of pressure levels per decade of pressure. This is sometimes a desirable option. The num_p and num_t arguments are the number of pressure levels to tabulate and the number of temperature levels to tabulate.

For some cryogenic cavitation problems, the requested (P,T) rectangle can cross onto the wrong side of the saturation boundary. A common example is a vapor table with T clipped low enough that no stable pure-vapor state exists above some pressure. In that situation the default behavior is to stop with an error because the requested phase does not exist in the specified temperature range.

If a rectangular table is still needed for cavitation_nist mixture-property lookups, the optional flag --extend_to_phase_boundary can be used. When enabled, the tool writes a single continuation row at pressure levels where the requested phase has no valid temperature window:

• For subcritical pressures, the continuation row is taken from the saturation boundary for that phase.

• If a later pressure level still needs continuation after the saturation boundary has terminated at the critical point, the tool reuses the last available continuation row as a fallback.

This option is intended for solver compatibility rather than for creating new thermodynamically exact pure-phase states. It should therefore be used deliberately and documented in the case setup.

An example command for a hydrogen vapor table that needs this continuation behavior is shown below.

<StreamInstallDir>/bin/cav_tab_tool liquidVapor --fluid_name HYDROGEN --t_min 16 --t_max 24 --p_min 20000 --p_max 6.0e6 --log_p --num_p 180 --num_t 160 --adaptive_t --t_rel_err 5e-4 --t_integrated_err 2e-3 --t_min_points 16 --t_max_points 450 --extend_to_phase_boundary

Saturation Tabulation

To generate a tabulated saturation file, one should pass the saturation argument to the tabulation tool. The presence of this argument enables a required set of arguments which are described below:

Saturation Tabulation Arguments

Argument	Description
t_min	The minimum temperature range value
t_max	The maximum temperature range value
num_t	The number of temperature levels
adaptive_t	Optional flag to enable adaptive temperature spacing
t_rel_err	Relative error scale for adaptive spacing
t_abs_err	Absolute error scale for adaptive spacing
t_integrated_err | Target normalized integrated error used to stop adaptive refinement
t_min_points	Minimum number of temperature points used by adaptive spacing
t_max_points	Maximum number of temperature points used by adaptive spacing

An example command to generate a saturation tabulation file for nitrogen in the temperature range of 65K-300K with 150 tabulation points is shown below.

<StreamInstallDir>/bin/cav_tab_tool saturation --fluid_name NITROGEN --t_min 65 --t_max 300 --num_t 150

In the example above, the saturation argument activates the saturation curve tabulation mode of the script. This mode has fewer arguments than the liquidVapor mode. The temperature range and number of points to include are the only arguments that are necessary. These arguments have the same meanings as they do for the liquidVapor mode of the tabulation tool that was discussed earlier.

Adaptive Temperature Spacing

When --adaptive_t is enabled, the tool refines the temperature grid by repeatedly splitting intervals that contribute the most interpolation error. The interpolation model remains piecewise linear; only the placement of points changes.

The adaptive controls have the following meanings:

t_rel_err

Relative error scale. A larger value is more permissive and usually produces fewer points.

t_abs_err

Absolute error scale. This is useful when a property crosses zero or is very small in magnitude.

t_integrated_err

Dimensionless global error target used to stop refinement. Smaller values add more points; larger values reduce points.

t_min_points

Minimum number of temperature points (safety floor).

t_max_points

Maximum number of temperature points (hard cap).

The integrated error used by the tool is a normalized quantity:

1. On each interval, evaluate the true property value at the midpoint.

2. Compare that midpoint value against linear interpolation from interval endpoints.

3. Normalize the midpoint error by t_abs_err + t_rel_err*scale, where scale is based on local property magnitude.

4. Use the largest normalized error across tabulated outputs for that interval.

5. Weight by interval length and sum over the full temperature range.

This produces a single dimensionless metric (the “integrated error”) that reflects both local interpolation quality and how much of the range is affected.

Practical tuning guidance:

• Start with --t_rel_err 0.02 --t_integrated_err 0.05 --t_min_points 6 --t_max_points 250.

• If files are too large, increase t_rel_err and/or t_integrated_err, or lower t_max_points.

• If interpolation quality is insufficient, decrease t_rel_err and/or t_integrated_err, or raise t_max_points.

• If t_max_points is reached, the tool reports a warning and stops refinement at the cap.

For current cavitation_nist code, saturation interpolation supports non-uniform temperature spacing, so adaptive spacing can be used for .sat files as well as .liq/.vap.

Using Cavitation Tabulation Tool

One must have version 11 of the REFPROP software to use the cavitation tabulation tool. A short tutorial on how to install the dependencies for the cavitation tabulation tool is presented below.

1. Download the REFPROP 11 software using the NIST downloader. This must be done on a Windows machine. The software will be downloaded into a directory named REFPROP.

2. Copy the REFPROP directory over to a Linux machine. The software can be put in a directory such as /home/<User>/software/refprop, where <User> would be your username on the machine). Note: do not rename the REFPROP source code directory because that may cause issues. For example, the path should look as follows: /home/<User>/software/refprop/REFPROP

3. Compile a local version of REFPROP on your Linux machine

   1. Do a recursive clone of the following repository: git clone --recursive https://github.com/usnistgov/REFPROP-cmake.git

   2. Go to the root of the cloned repository: cd REFPROP-cmake

   3. Create a build directory: mkdir build

   4. Change to the build directory: cd build

   5. Run cmake to configure the build: cmake .. -DCMAKE_BUILD_TYPE=Release -DREFPROP_FORTRAN_PATH=/home/<User>/software/refprop/REFPROP/FORTRAN

   6. Build the project: cmake --build .

4. Once the build is complete, a librefprop.so file should be found in the build directory. This shared library contains all the REFPROP functions that the Python wrapper will call. Move this file to the REFPROP directory created in step 2.

   1. Move the shared library file: mv librefprop.so /home/<User>/software/refprop/REFPROP

5. Set the environment variable RPPREFIX to point to the REFPROP directory that contains the shared library file.

   1. Add the following line to your .bashrc file: export RPPREFIX=/home/<User>/software/refprop/REFPROP.

   2. Make sure to source your .bashrc file after adding this line: source ~/.bashrc.

6. Clone the REFPROP wrappers directory into the /home/<User>/software/refprop directory using the following command: git clone https://github.com/usnistgov/REFPROP-wrappers

7. Create a Python3 environment with the REFPROP interface installed. This will provide a portable and dependable Python environment for running the tabulation tool. Go to your refprop directory and follow the steps below:

   1. Create a virtual environment: python3 -m venv refprop-env

   2. Activate the environment: source refprop-env/bin/activate

   3. Install numpy: pip install numpy

   4. Navigate to the REFPROP wrappers directory: cd /home/<User>/software/refprop/REFPROP-wrappers/wrappers/python/ctypes

   5. Install the python wrapper into the environment: python setup.py install

   6. To leave the virtual environment type: deactivate

One you are in the python environment, you can run the tabulation script and it should execute using calls to the REFPROP library.

References

[MeFB1998]

3. 12. Merkle, J. Z. 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.

[HoAU2007]

1. Hosangadi, V. Ahuja and R. J. Ungewitter, “Analysis of Thermal Effects in Cavitating Liquid Hydrogen Inducers,” Journal of Propulsion and Power, vol. 23, no. 6, pp. 1225-1234, 2007.

[GeZB2004]

1. 7. Gerber, P. J. Zwart and T. Belamri, “A Two-Phase Flow Model for Predicting Cavitation Dynamics,” International Conference on Multiphase Flow, Yokohama, Japan, 2004.

[ScSS2008]

7. 8. Schnerr, I. H. Sezal and S. J. Schmidt, “Numerical investigation of three-dimensional cloud cavitation with special emphasis on collapse induced shock dynamics,” Physics of Fluids, vol. 20, 2008.