Simulation of Cavitating Flows

The cavitation_nist module models liquid, vapor, and non-condensable gas (NCG) at a shared pressure, temperature, and velocity. Liquid and vapor properties come from REFPROP tables. See Homogeneous Cavitation Modeling in the Theory Guide for the formulation.

Basic Setup

Load cavitation_nist, select compressible flow and the temperature energy equation, and supply the property tables and source-model parameters. Total-enthalpy and total-energy forms are not supported.

These entries configure cavitation for a laminar nitrogen case. Add the mesh and time-stepping controls for your case, and match the boundary names to the mesh. Values in this chapter’s examples illustrate syntax; defaults are listed separately.

loadModule: cavitation_nist
flowRegime: laminar
flowCompressibility: compressible
energyEquationOptions: <form=temperature>

liquid_model: NITROGEN
vapor_model: NITROGEN
transport_model: module

cavitationEquationOptions: <source=Zwart>
ZwartSourceParameters: <RB=1e-6, aNuc=5e-4, fVap=50, fCond=0.01>

initialCondition: <p=200000 Pa, T=80 K, v=[10 m/s,0,0],
                   vapor_y=0, NCG_y=0>
boundary_conditions: <
    Inlet=subsonicInlet(v=10 m/s, T=80 K, vapor_y=0, NCG_y=0),
    Outlet=fixedPressureOutlet(p=200000 Pa),
    Wall=noslip
>
plot_output: vapor_y,vapor_vof,NCG_y

liquid_model and vapor_model are filename prefixes. This example requires NITROGEN.liq, NITROGEN.vap, and NITROGEN.sat; the saturation file uses the liquid prefix. transport_model: module selects viscosity and conductivity from the tables. For table generation, see Property-Table Generation.

vapor_y and NCG_y are mass fractions, not volume fractions. Both must be nonnegative and their sum must not exceed one. Use NCG_y=0 when no non-condensable gas is present.

Non-Condensable Gas

NCG properties default to air. To change them, specify nonCondensableGasProperties; values must be positive and omitted entries retain their defaults:

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

NCG properties

Entry	Meaning and units	Default
m	Molar mass, g/mol	28.9647
cp	Specific heat at constant pressure, J/(kg K)	1005
mu	Dynamic viscosity, kg/(m s)	1.805e-5
kcond	Thermal conductivity, W/(m K)	0.02476

Initial and Boundary Conditions

Specify vapor_y and NCG_y with the initial flow state and at subsonicInlet, supersonicInlet, and inflow boundaries. For spatially varying initialization, see Appendix: Initial Conditions.

Open boundaries

Boundary	Composition and incoming flow
fixedPressureOutlet	Extrapolates composition, including during backflow. Use only where the normal flow remains subsonic at the local mixture sound speed.
outflow	Extrapolates composition and blocks inward flow. Does not prescribe an incoming reservoir state.
farfield	Requires vapor_y; NCG_y defaults to zero. Uses the reference state for incoming flow and the interior state for outgoing flow.

For a boundary that admits incoming reservoir fluid, use farfield:

Open=farfield(p=200000 Pa, T=80 K, v=[10 m/s,0,0], vapor_y=0, NCG_y=0)

For subsonic outflow and farfield, boundary temperature is recovered from the mixture EOS within TclipMin and TclipMax. Failure to find a converged temperature stops the run. A farfield reference temperature need not equal the resulting face temperature. Choose bounds consistent with the property tables; widening them does not extend the tables’ physical validity. See Appendix: Boundary Conditions for the remaining boundary inputs.

Cavitation Models

Select Merkle, Zwart, or SauerSchnerr through source in cavitationEquationOptions.

Merkle

Merkle rates depend on pressure departure from vapor pressure and the reference scales below. All five entries are required.

cavitationEquationOptions: <source=Merkle>
MerkleSourceParameters: <tauVap=1, tauCond=80, L=0.7, V=1.5, rho=800>

MerkleSourceParameters

Entry	Meaning	Units
tauVap	Vaporization time-scale parameter; larger values weaken vaporization	Dimensionless
tauCond	Condensation time-scale parameter; larger values weaken condensation	Dimensionless
L	Reference length	m
V	Reference speed	m/s
rho	Reference density	kg/m3

Zwart

Zwart uses a prescribed bubble radius and nuclei volume fraction. All four entries are required; the basic setup shows their syntax.

ZwartSourceParameters

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

Sauer-Schnerr

Sauer-Schnerr derives bubble size from vapor content and nuclei density. For a fixed nuclei density:

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

SauerSchnerrSourceParameters

Entry	Meaning and units	Default
rNuc	Seed radius for vaporization, m; finite and positive	1e-5
n	Fixed nuclei number density, m-3; finite and positive	Required

To transport nuclei instead of prescribing n, see Transported Nuclei.

Source Scaling

An optional multiplier modifies the source below cavitationMinPressure. The multiplier is one at that pressure and reaches cavitationMaxSourceFactor at zero pressure. Factors above one amplify phase change; factors between zero and one damp it.

Source scaling inputs

Input	Meaning	Default
cavitationMaxSourceFactor	Dimensionless factor at zero pressure; 1 disables scaling	1
cavitationSourceFactorExponent	Dimensionless curve exponent; nonzero. Larger positive values retain the zero-pressure factor over more of the pressure interval.	1
cavitationMinPressure	Activation pressure in Pa; nonzero. This is not a pressure floor.	1

Numerical Options

Add solver settings to the selected source in cavitationEquationOptions:

cavitationEquationOptions: <source=Zwart, linearSolver=SGS,
                            relaxationFactor=0.5, maxIterations=5>

Cavitation equation settings

Entry	Meaning	Default
linearSolver	SGS, FSGS, LSGS, PETSC, or HYPRE; see Appendix: Linear Solvers	SGS
relaxationFactor	Equation relaxation, greater than zero and at most one	0.5
maxIterations	Maximum linear iterations	5
conservativeFluxRescale	on limits higher-order phase fluxes; see below	off

The following are top-level run control inputs:

Transport and pressure settings

Input	Meaning	Default
cavitation	on or off; off disables vapor sources, retaining phase transport	on
cavitationInviscidFlux	Vapor and NCG transport: FOU or SOU	SOU
AlphaViscosity	Artificial diffusion coefficient for vapor and NCG, kg/(m s)	1e-5
liquid_rho, vapor_rho, NCG_rho	Reference densities for pressure-equation normalization, kg/m3. With liquid_rho=-1, local phase densities are used. To prescribe references, set all three to positive values.	-1

Conservative Flux Rescaling

SOU bounds reconstructed face fractions. conservativeFluxRescale=on also limits the higher-order flux correction using available phase mass. This can reduce cell clipping at the cost of additional numerical diffusion. The selected reconstruction limiter still applies.

For the Zwart setup above:

cavitationEquationOptions: <source=Zwart, conservativeFluxRescale=on>
cavitationInviscidFlux: SOU
timeIntegrator: BDF
PDEContinuityAdded: yes
AlphaViscosity: 0
explicitCorrectorNewtonIteration: off

Use a fixed mesh, BDF or BDF2, and the settings shown. All three source models are supported, including transported nuclei and cavitation: off. SIMPLE/SIMPLEC and PIMPLE are supported; moving meshes and the Newton explicit corrector are not.

BDF2 requires a nonnegative weighted mass history for vapor, NCG, and liquid. An invalid history stops the run; reduce the timestep or restart with BDF. PIMPLE also stops if its explicit source update exceeds available phase mass; reduce the timestep. Neither failure triggers an automatic retry.

Convergence

Both phase equations must converge, including with cavitation: off. Use vapor and NCG entries in the convergence tolerances:

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

The absolute phase residual is the sum of absolute cell imbalances in kg/s. The relative residual is measured against the first positive residual of that phase in the current timestep. Either tolerance can satisfy convergence, with an allowance for roundoff. Entries must be finite and nonnegative. Omitted entries inherit default, or 1e-30 if it is absent. An active legacy convergenceTolerance overrides the absolute tolerance and sets the relative tolerance to 1e-30.

An unconverged final phase state stops the run, including at maxIterationsPerTimeStep. Stream does not accept that timestep or reduce its size automatically. Check residuals, phase bounds, timestep, and spatial discretization before restarting.

Turbulence Corrections

DFR Eddy-Viscosity Correction

DFR reduces turbulent eddy viscosity in liquid-vapor mixtures. To add it to the basic setup, use turbulent flow and SST 2003 RANS:

flowRegime: turbulent
turbulenceEquationOptions: <model=menterSST2003>
cavitationEquationOptions: <source=Zwart, turbulenceCorrection=DFR>
DFRTurbulenceCorrectionParameters: <n=10>

Supply turbulence initial and boundary data as described in Simulation of Turbulent Flows. The exponent n defaults to 10 and must be finite and at least one; n=1 leaves eddy viscosity unmodified. Omit turbulenceCorrection to disable DFR. Other turbulence models and des, DESLimiter, and rhok are unsupported with DFR.

DFR supports all cavitation source models and cavitation: off. It leaves the mixture EOS unchanged. Prescribed boundary eddy viscosity, including muRatio, remains an override. Its exponent requires calibration for the flow of interest.

Vapor-Pressure Correction

turbulentVapourPressureCorrection: on raises the phase-change pressure threshold to account for turbulent pressure fluctuations. The default is off. This option is independent of DFR and applies to turbulent flow.

Transported Nuclei

For Sauer-Schnerr, replace fixed n with a transported nuclei density using nucleationTransport. nb is nuclei per cubic metre of mixture; it is not interchangeable with the fixed model’s reference-volume n.

cavitationEquationOptions: <source=SauerSchnerr>
SauerSchnerrSourceParameters: <rNuc=3e-6, nucleationTransport>
nucleationEquationOptions: <linearSolver=SGS, relaxationFactor=0.5,
                            maxIterations=5>
nucleationInviscidFlux: FOU
timeIntegrator: BDF

Supply finite, nonnegative nb in the initial state and at every inlet and farfield boundary. For example, extend the basic setup with:

initialCondition: <p=200000 Pa, T=80 K, v=[10 m/s,0,0],
                   vapor_y=0, NCG_y=0, nb=1e8>
boundary_conditions: <
    Inlet=subsonicInlet(v=10 m/s, T=80 K, vapor_y=0, NCG_y=0, nb=1e8),
    Outlet=fixedPressureOutlet(p=200000 Pa),
    Wall=noslip
>

Do not specify fixed n as well. Outlets extrapolate nb; an optional outlet nb supplies a backflow reservoir. Walls have zero advective number flux. Only BDF and BDF2 are supported. Start with FOU for nuclei.

Nuclei controls

Input	Meaning	Default
nucleationInviscidFlux	Top-level input: FOU or SOU	SOU
nbViscosity	Top-level artificial diffusivity, m2/s; nonnegative	0
absoluteTolerance	In nucleationEquationOptions: absolute number residual, nuclei/s	1e-3
relativeTolerance	In nucleationEquationOptions: relative residual; at least zero and less than one	1e-8

Solver controls default to linearSolver=SGS, relaxationFactor=0.5, and maxIterations=5. Both tolerances must be finite and nonnegative, and at least one must be positive. A materially negative nb or an unconverged number equation at the iteration limit stops the run. Reduce the timestep or use a more bounded discretization when bounds fail.

Nuclei restarts require nbMixture_hdf5 and, for BDF2, nbMixturePrev_hdf5. Restarts containing only the older nb_hdf5 format cannot be used; reinitialize nuclei.

Bubble Growth

Select bubbleGrowth in nucleationEquationOptions:

Bubble growth choices

Value	Effect
rayleigh (default)	Inertial growth and collapse from pressure departure
reducedRP	Adds surface tension, liquid viscosity, and isothermal gas pressure; requires positive liquid viscosity and nonnegative surface tension

rNuc sets a minimum radius for the rate calculation; it adds no phase mass. Zero nb gives zero phase-change rate. reducedRP is experimental: it does not resolve radial inertia, collapse rebound, or acoustic emission. Validate it for the intended application.

Population Sources

Select population in nucleationEquationOptions. Breakup and coalescence change bubble count and representative size without adding phase mass.

Population choices and parameters

Value	Parameters in nucleationEquationOptions
none (default)	No breakup or coalescence
constant	breakupRate in s-1 and coalescenceKernel in m3/s; both default to zero and must be nonnegative
turbulentSingleSize	breakupMultiplier and coalescenceMultiplier default to one; nonnegative, with zero disabling that contribution. filmThicknessRatio defaults to 1e4 and must exceed one.

turbulentSingleSize requires menterSST2003 RANS without des, DESLimiter, or rhok. It models one representative bubble size, not a size distribution. Its rates require validation for the liquid and flow.

Classical Nucleation

nucleation=classical in nucleationEquationOptions adds homogeneous and wall nucleation. The default is none. These pressure- and wetting-sensitive laws do not replace measured inlet nuclei.

Classical nucleation inputs

Input	Meaning	Default
SpeciesMolecularMass	In nucleationEquationOptions: molar mass in kg/kmol; positive	Required
thetaB	Top-level bubble contact angle in degrees; greater than zero and at most 180	170
nbSourceClipMax	Top-level birth-rate limit, nuclei/(m3 s); positive	1e30

Output and Diagnostics

Add the following fields to plot_output as needed. All are optional; units below describe the cell values before nodal interpolation.

Mixture and phase fields

Variable	Meaning	Units
liquid_y, vapor_y, NCG_y	Mass fractions	Dimensionless
liquid_vof, vapor_vof, NCG_vof	Volume fractions	Dimensionless
rhoL, rhoV	Saturation liquid and vapor densities	kg/m3
pV	Vapor pressure used by the source	Pa
bt	Bulk expansion coefficient times temperature	Dimensionless
F	Source scaling factor	Dimensionless
rhoFactor	DFR eddy-viscosity multiplier	Dimensionless
vapor_yResidual, NCG_yResidual	Cell equation imbalances	kg/s
phaseFluxReduction	Higher-order correction reduction: zero is unrestricted, one removes it	Dimensionless
AlphaProductionRate	Cell-integrated vapor source; positive for evaporation, negative for condensation	kg/s
AlphaProductionTimeScale	Available phase mass divided by source magnitude; capped at 1e30 for zero or very small rates	s

AlphaProductionTimeScale is a depletion-time diagnostic, not a stable timestep estimate. It and AlphaProductionRate use maximum-magnitude cell-to-node interpolation. phaseFluxReduction uses the largest adjacent cell value and is zero when rescaling is off.

Nuclei fields

Variable	Meaning and availability	Units
nb	Nuclei density; transported or prescribed	m-3
nucleiRadius	Representative radius; transported nuclei	m
nucleiGasPressure	Bubble gas pressure; transported nuclei	Pa
nucleiBreakupRate	Breakup frequency; transported nuclei	s-1
nucleiCoalescenceKernel	Coalescence kernel; transported nuclei	m3/s
nbResidual	Cell number-equation imbalance; transported nuclei	nuclei/s
Jhom	Homogeneous nucleation rate; classical nucleation	nuclei/(m3 s)
Jsrc	Cell-integrated nucleation rate before nbSourceClipMax; classical nucleation	nuclei/s

Flux Diagnostics

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

At print_freq, the CFLUX: line reports cycle, iteration, cell count, limited-cell count, mean correction reduction, materially clipped-cell count, and absolute vapor and NCG clipping adjustments in kg. These are iteration diagnostics; summing adjustments over iterations does not measure physical mass loss.

The energy column in flux_<boundary>.dat includes convected internal and kinetic energy relative to the grid, plus heat transfer. It excludes pressure and viscous work, so it is not a complete total-energy balance. Liquid/vapor energy datums come from the tables; NCG internal energy is zero at zero kelvin.

The component mass diagnostic reports the signed mass-balance defect in kg/s using predictor fluxes and the configured time history. It should approach zero as continuity converges; it does not replace an accepted-state mass budget.

Property-Table Generation

cav_tab_tool, installed in Stream’s bin/ directory, generates the .liq, .vap, and .sat files for cavitation_nist. The separate tabulated_data utility generates EoSTables/ data for other EOS models.

Prerequisites

Install REFPROP and its shared library using the REFPROP build instructions. Set RPPREFIX to the REFPROP directory containing the library and fluid data. Install NumPy and the Python wrapper in the Python environment used to run the tool:

python3 -m venv refprop-env
. refprop-env/bin/activate
python -m pip install numpy ctREFPROP
export RPPREFIX=/path/to/REFPROP

Generate Tables

Run both commands to create NITROGEN.liq, NITROGEN.vap, and NITROGEN.sat in the current directory. Choose ranges that cover the case:

cav_tab_tool liquidVapor --fluid_name NITROGEN \
    --t_min 65 --t_max 130 --p_min 50000 --p_max 500000 \
    --num_t 200 --num_p 200
cav_tab_tool saturation --fluid_name NITROGEN \
    --t_min 65 --t_max 125 --num_t 150

Tabulation arguments

Argument	Meaning	Default
--fluid_name	REFPROP fluid name; must match the installed fluid database	Required
--t_min, --t_max	Temperature bounds in K; minimum less than maximum	Required
--num_t	Uniform temperature points; at least two	200
--p_min, --p_max	Pressure bounds in Pa; minimum less than maximum. liquidVapor only.	Required
--num_p	Pressure points; at least two. liquidVapor only.	200
--log_p	Uniform spacing in log pressure; requires positive bounds. liquidVapor only.	Off

Temperature windows follow the stable phase at each pressure. The liquid window can extend below --t_min to the EOS minimum. Saturation tables stop below the critical temperature even if --t_max is higher.

Adaptive Spacing

Add --adaptive_t to either command to refine temperature spacing where properties are curved. Interpolation remains piecewise linear. These controls replace the uniform --num_t count:

Adaptive temperature controls

Argument	Effect	Default
--t_rel_err	Relative error scale; smaller values generally add points	1e-3
--t_abs_err	Absolute error scale for properties near zero	0
--t_integrated_err	Dimensionless integrated error target; smaller values add points	1e-2
--t_min_points	Minimum temperature points; at least two	4
--t_max_points	Maximum temperature points; at least --t_min_points	2000

Error scales must be nonnegative, with at least one positive. The integrated error target must be positive. Reaching the point cap before meeting the target produces a warning; increase the cap or relax the error controls.

Phase-Boundary Continuation

If a pressure level has no valid temperature window for a phase, the tool stops by default. Add --extend_to_phase_boundary to liquidVapor to write a single saturation-boundary row at that pressure instead. Beyond the critical point, continuation reuses the last available boundary row; without one it fails.

Continuation supplies numerical lookup values, not valid pure-phase states outside the phase boundary. Record its use in the case setup.