pTg phase-change activation is blocked by the vapor-mass threshold in the six-equation cavitation model

Summary

I am using the MFC six-equation model (model_eqns = 3) with relax = 'T' and relax_model = 6 (pTg equilibrium) to simulate a high-speed submerged waterjet and the associated cavitation cloud.

I believe there is a potential issue with the current pTg phase-change activation criterion in src/common/m_phase_change.fpp.

The current criterion requires both the liquid and vapor partial masses to satisfy:

m_l > mixM*rM
m_v > mixM*rM

before s_infinite_ptg_relaxation_k() is called.

With the default:

mixM = 1.0e-8_wp

this prevents a physically small initial vapor seed from triggering liquid-to-vapor phase change, even when the liquid state has entered a cavitation regime.

Relevant source logic

In the current source, the pTg relaxation is guarded by a condition of the form:

if ((relax_model == 6) .and. &
    ((m_l > mixM*rM) .and. (m_v > mixM*rM)) .and. &
    (pS < pCr) .and. (TS < TCr)) then

    call s_infinite_ptg_relaxation_k(...)
end if

Therefore, for a trace-vapor cell, pTg relaxation cannot start unless the vapor partial mass has already exceeded mixM*rM.

This appears problematic for cavitation initiation because the vapor mass may need to grow from a very small vapor seed in the first place.

The pTg Newton solver itself appears to conserve the total reacting mass,

mT = m_l + m_v

and solve for m_l, with

m_v = mT - m_l

so the m_v > mixM*rM requirement appears to be an activation/safeguard criterion outside the pTg solver rather than a mathematical requirement of the pTg equilibrium solve itself.

My test case

The case is a submerged high-speed medical waterjet:

• MFC version: v5.7.0-37-gfef3c276
• model_eqns = 3 (six-equation model)
• relax = 'T'
• relax_model = 6
• liquid initial density: approximately 998 kg/m3
• initial vapor volume fraction: alpha_v = 1e-8
• initial vapor density: approximately 0.0173 kg/m3
• saturation pressure at the operating temperature: approximately 2340 Pa
• tank pressure: approximately 1e5 Pa
• jet velocity scale: approximately 300–320 m/s

With:

mixM = 1e-8

the activation threshold is approximately:

mixM*rM ≈ 1e-8 × 998 ≈ 9.98e-6 kg/m3

However, the initial vapor partial mass is only:

m_v = alpha_v * rho_v
    ≈ 1e-8 × 0.0173
    ≈ 1.73e-10 kg/m3

Therefore the initial vapor partial mass is approximately 60,000 times below the current activation threshold.

Observed cell

I tracked a cell that eventually develops a NaN state:

cell = (49, 246, 61)
position ≈ (-0.0691, -0.9272, 0.0349) mm
radial position ≈ 77.4 µm

This location is essentially at the nozzle lip / shear-layer region, where cavitation inception is expected in this configuration.

The recorded conservative variables for this cell are:

Frame 23:
a1r1 = 9.549e+02
a2r2 = 1.654e-10
al1  = 9.580e-01
al2  = 4.202e-02

Frame 24:
a1r1 = 9.529e+02
a2r2 = 1.650e-10
al1  = 9.561e-01
al2  = 4.394e-02

Frame 25:
a1r1 = 9.510e+02
a2r2 = 1.647e-10
al1  = 9.541e-01
al2  = 4.588e-02

The important observation is:

vapor partial mass:
1.654e-10 → 1.650e-10 → 1.647e-10 kg/m3

which remains essentially constant.

At the same time, the vapor volume fraction increases:

alpha_v:
0.04202 → 0.04394 → 0.04588

This corresponds to a drastic decrease in reconstructed vapor density:

rho_v ≈ 3.94e-9 → 3.75e-9 → 3.59e-9 kg/m3

The volume fractions themselves remain approximately normalized:

alpha_l + alpha_v ≈ 1

Therefore, this does not appear to be primarily a volume-fraction closure problem.

Instead, the behavior is consistent with a trace vapor seed becoming extremely rarefied while its vapor mass does not significantly increase.

What happens when I reduce mixM

As a diagnostic experiment, I reduced:

mixM = 1e-8

to:

mixM = 1e-14

This allowed many trace-vapor cells to satisfy the pTg activation condition.

However, the pre-processing stage then became extremely slow / appeared to hang after reporting:

initial condition might have been altered due to enforcement of pTg-equilibrium (relax = "T" activated)

This suggests that globally reducing mixM is not a satisfactory solution because it causes a very large number of nearly pure-liquid cells containing tiny vapor seeds to become pTg candidates.

I therefore do not want to solve the problem simply by changing mixM.

Question / potential issue

Could the pTg activation criterion be decoupled from the mixM threshold for cavitation initiation?

In particular, is the current requirement

m_v > mixM*rM

intended to prevent pTg from initiating phase change from a physically small vapor seed?

If so, how is cavitation nucleation/growth from a trace vapor seed intended to occur in the six-equation model when:

m_v << mixM*rM

but the liquid state has already entered a thermodynamically unstable or metastable cavitation regime?

One possible approach would be to retain mixM as a numerical safeguard for near-single-phase states, while providing a separate cavitation activation criterion for trace-vapor states, for example a criterion based on the liquid thermodynamic state such as:

p_l < p_sat(T_l)

followed by the existing pTg equilibrium solver.

I am not suggesting that this particular criterion must be the final implementation; rather, I would like to clarify whether MFC currently has an intended mechanism for this situation.

Expected behavior

For a cell containing a very small vapor seed:

m_v << mixM*rM

when the liquid enters a cavitation regime, I would expect the phase-change model to be able to initiate/grow the vapor phase rather than requiring the vapor phase to first reach the mass threshold that the phase-change mechanism itself is supposed to help produce.

At the same time, I would expect the existing mixM safeguard to remain available so that pTg is not unnecessarily executed in every trace-vapor cell throughout the domain.

Request

Could the developers please clarify:

1. Is the m_v > mixM*rM condition intentionally required for pTg activation in the six-equation cavitation model?
2. If intentional, what is the intended mechanism for cavitation initiation from a trace vapor seed with m_v << mixM*rM?
3. Is there an existing recommended configuration or source-level mechanism for this situation?
4. If there is no such mechanism, would it be appropriate to introduce a separate cavitation/nucleation activation criterion that does not require the vapor mass to exceed mixM*rM, while retaining mixM for numerical protection?

I can provide the complete case and additional cell/time-history data if needed.