From 9567bc0d4b13762917eb96a0fc2a8251d248d06c Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Tue, 15 Nov 2022 12:52:41 +0000
Subject: [PATCH] cavitationModels: Made insensitive to phase order

The cavitation models used by the interFoam solver and the
compressibleVoF solver module can now be applied regardless of the
ordering of the liquid and vapour phases. A "liquid" keyword is now
required in the model specification in order to control which phase is
considered to be the condensed liquid state. Previously the liquid phase
was assumed to be the first of the two phases.
---
 .../fvModels/VoFCavitation/VoFCavitation.C    | 24 ++++----
 .../fvModels/VoFCavitation/VoFCavitation.C    | 24 ++++----
 .../compressibleCavitation/Kunz/Kunz.C        | 30 +++++-----
 .../compressibleCavitation/Kunz/Kunz.H        |  8 +--
 .../compressibleCavitation/Merkle/Merkle.C    | 19 +++---
 .../compressibleCavitation/Merkle/Merkle.H    |  8 +--
 .../SchnerrSauer/SchnerrSauer.C               | 40 ++++++-------
 .../SchnerrSauer/SchnerrSauer.H               | 12 ++--
 .../cavitationModel/cavitationModel.C         |  1 +
 .../cavitationModel/cavitationModel.H         | 58 ++++++++++++++++---
 .../compressibleTwoPhases.H                   | 12 ++++
 .../incompressibleCavitation/Kunz/Kunz.C      | 30 +++++-----
 .../incompressibleCavitation/Kunz/Kunz.H      |  8 +--
 .../incompressibleCavitation/Merkle/Merkle.C  | 18 +++---
 .../incompressibleCavitation/Merkle/Merkle.H  |  8 +--
 .../SchnerrSauer/SchnerrSauer.C               | 40 ++++++-------
 .../SchnerrSauer/SchnerrSauer.H               | 12 ++--
 .../cavitationModel/cavitationModel.C         |  1 +
 .../cavitationModel/cavitationModel.H         | 56 ++++++++++++++++--
 .../incompressibleTwoPhases.H                 | 11 +++-
 .../twoPhaseMixture/twoPhases.H               | 31 +++++++++-
 .../ballValve/0/{alpha.water => alpha.vapour} |  4 +-
 .../ballValve/constant/fvModels               |  2 +
 .../ballValve/constant/phaseProperties        |  2 +-
 .../ballValve/system/fvSolution               |  2 +-
 .../interFoam/RAS/propeller/constant/fvModels |  6 ++
 .../cavitatingBullet/constant/fvModels        |  6 ++
 27 files changed, 307 insertions(+), 166 deletions(-)
 rename tutorials/modules/compressibleVoF/ballValve/0/{alpha.water => alpha.vapour} (95%)

diff --git a/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C b/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
index f2ac15c4bd..a5bc25f4f6 100644
--- a/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
+++ b/applications/solvers/modules/fluid/compressibleVoF/fvModels/VoFCavitation/VoFCavitation.C
@@ -96,21 +96,21 @@ void Foam::fv::compressible::VoFCavitation::addSup
 
     if (fieldName == alphaName_)
     {
-        volScalarField::Internal alphalCoeff
+        const volScalarField::Internal alpha1Coeff
         (
             1.0/mixture_.rho1()()
           - mixture_.alpha1()()*(1.0/mixture_.rho1()() - 1.0/mixture_.rho2()())
         );
 
-        const Pair<tmp<volScalarField::Internal>> mDotAlphal
+        const Pair<tmp<volScalarField::Internal>> mDot12Alpha
         (
-            cavitation_->mDotAlphal()
+            cavitation_->mDot12Alpha()
         );
 
-        const volScalarField::Internal vDotcAlphal(alphalCoeff*mDotAlphal[0]());
-        const volScalarField::Internal vDotvAlphal(alphalCoeff*mDotAlphal[1]());
+        const volScalarField::Internal vDot1Alpha(alpha1Coeff*mDot12Alpha[0]());
+        const volScalarField::Internal vDot2Alpha(alpha1Coeff*mDot12Alpha[1]());
 
-        eqn += fvm::Sp(vDotvAlphal - vDotcAlphal, eqn.psi()) + vDotcAlphal;
+        eqn += fvm::Sp(-vDot2Alpha - vDot1Alpha, eqn.psi()) + vDot1Alpha;
     }
 }
 
@@ -140,17 +140,17 @@ void Foam::fv::compressible::VoFCavitation::addSup
             1.0/mixture_.rho1()() - 1.0/mixture_.rho2()()
         );
 
-        const Pair<tmp<volScalarField::Internal>> mDotP
+        const Pair<tmp<volScalarField::Internal>> mDot12P
         (
-            cavitation_->mDotP()
+            cavitation_->mDot12P()
         );
 
-        const volScalarField::Internal vDotcP(pCoeff*mDotP[0]);
-        const volScalarField::Internal vDotvP(pCoeff*mDotP[1]);
+        const volScalarField::Internal vDot1P(pCoeff*mDot12P[0]);
+        const volScalarField::Internal vDot2P(pCoeff*mDot12P[1]);
 
         eqn +=
-            (vDotvP - vDotcP)*(cavitation_->pSat() - rho*gh)
-          - fvm::Sp(vDotvP - vDotcP, eqn.psi());
+            (vDot2P - vDot1P)*(cavitation_->pSat() - rho*gh)
+          - fvm::Sp(vDot2P - vDot1P, eqn.psi());
     }
 }
 
diff --git a/applications/solvers/multiphase/interFoam/fvModels/VoFCavitation/VoFCavitation.C b/applications/solvers/multiphase/interFoam/fvModels/VoFCavitation/VoFCavitation.C
index 7723e4362c..cdf68f8d0c 100644
--- a/applications/solvers/multiphase/interFoam/fvModels/VoFCavitation/VoFCavitation.C
+++ b/applications/solvers/multiphase/interFoam/fvModels/VoFCavitation/VoFCavitation.C
@@ -95,21 +95,21 @@ void Foam::fv::VoFCavitation::addSup
 
     if (fieldName == alphaName_)
     {
-        volScalarField::Internal alphalCoeff
+        const volScalarField::Internal alpha1Coeff
         (
             1.0/mixture_.rho1()
           - mixture_.alpha1()()*(1.0/mixture_.rho1() - 1.0/mixture_.rho2())
         );
 
-        const Pair<tmp<volScalarField::Internal>> mDotAlphal
+        const Pair<tmp<volScalarField::Internal>> mDot12Alpha
         (
-            cavitation_->mDotAlphal()
+            cavitation_->mDot12Alpha()
         );
 
-        const volScalarField::Internal vDotcAlphal(alphalCoeff*mDotAlphal[0]());
-        const volScalarField::Internal vDotvAlphal(alphalCoeff*mDotAlphal[1]());
+        const volScalarField::Internal vDot1Alpha(alpha1Coeff*mDot12Alpha[0]());
+        const volScalarField::Internal vDot2Alpha(alpha1Coeff*mDot12Alpha[1]());
 
-        eqn += fvm::Sp(vDotvAlphal - vDotcAlphal, eqn.psi()) + vDotcAlphal;
+        eqn += fvm::Sp(-vDot2Alpha - vDot1Alpha, eqn.psi()) + vDot1Alpha;
     }
 }
 
@@ -139,17 +139,17 @@ void Foam::fv::VoFCavitation::addSup
             1.0/mixture_.rho1() - 1.0/mixture_.rho2()
         );
 
-        const Pair<tmp<volScalarField::Internal>> mDotP
+        const Pair<tmp<volScalarField::Internal>> mDot12P
         (
-            cavitation_->mDotP()
+            cavitation_->mDot12P()
         );
 
-        const volScalarField::Internal vDotcP(pCoeff*mDotP[0]);
-        const volScalarField::Internal vDotvP(pCoeff*mDotP[1]);
+        const volScalarField::Internal vDot1P(pCoeff*mDot12P[0]);
+        const volScalarField::Internal vDot2P(pCoeff*mDot12P[1]);
 
         eqn +=
-            (vDotvP - vDotcP)*(cavitation_->pSat() - rho*gh)
-          - fvm::Sp(vDotvP - vDotcP, eqn.psi());
+            (vDot2P - vDot1P)*(cavitation_->pSat() - rho*gh)
+          - fvm::Sp(vDot2P - vDot1P, eqn.psi());
     }
 }
 
diff --git a/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.C b/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.C
index 11c8b13a71..bb81de1125 100644
--- a/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.C
+++ b/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.C
@@ -65,55 +65,53 @@ Foam::compressible::cavitationModels::Kunz::Kunz
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::compressible::cavitationModels::Kunz::mDotAlphal() const
+Foam::compressible::cavitationModels::Kunz::mDotcvAlphal() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
-    const volScalarField::Internal mcCoeff_(Cc_*phases_.rho2()()/tInf_);
+    const volScalarField::Internal mcCoeff_(Cc_*rhov()/tInf_);
     const volScalarField::Internal mvCoeff_
     (
-        Cv_*phases_.rho2()()/(0.5*phases_.rho1()()*sqr(UInf_)*tInf_)
+        Cv_*rhov()/(0.5*rhol()*sqr(UInf_)*tInf_)
     );
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        mcCoeff_*sqr(limitedAlpha1)
+        mcCoeff_*sqr(limitedAlphal)
        *max(p - pSat(), p0_)/max(p - pSat(), 0.01*pSat()),
-
-        mvCoeff_*min(p - pSat(), p0_)
+       -mvCoeff_*min(p - pSat(), p0_)
     );
 }
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::compressible::cavitationModels::Kunz::mDotP() const
+Foam::compressible::cavitationModels::Kunz::mDotcvP() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
-    const volScalarField::Internal mcCoeff_(Cc_*phases_.rho2()()/tInf_);
+    const volScalarField::Internal mcCoeff_(Cc_*rhov()/tInf_);
     const volScalarField::Internal mvCoeff_
     (
-        Cv_*phases_.rho2()()/(0.5*phases_.rho1()()*sqr(UInf_)*tInf_)
+        Cv_*rhov()/(0.5*rhol()*sqr(UInf_)*tInf_)
     );
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        mcCoeff_*sqr(limitedAlpha1)*(1.0 - limitedAlpha1)
+        mcCoeff_*sqr(limitedAlphal)*(1.0 - limitedAlphal)
        *pos0(p - pSat())/max(p - pSat(), 0.01*pSat()),
-
-        (-mvCoeff_)*limitedAlpha1*neg(p - pSat())
+        (-mvCoeff_)*limitedAlphal*neg(p - pSat())
     );
 }
 
diff --git a/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.H b/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.H
index 5525e709a3..6f08d09f16 100644
--- a/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.H
+++ b/src/twoPhaseModels/compressibleCavitation/Kunz/Kunz.H
@@ -102,13 +102,13 @@ public:
     // Member Functions
 
         //- Return the mass condensation and vaporisation rates as a
-        //  coefficient to multiply (1 - alphal) for the condensation rate
-        //  and a coefficient to multiply  alphal for the vaporisation rate
-        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const;
+        //  coefficient to multiply alphav for the condensation rate and a
+        //  coefficient to multiply alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvAlphal() const;
 
         //- Return the mass condensation and vaporisation rates as coefficients
         //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField::Internal>> mDotP() const;
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvP() const;
 
         //- Correct the Kunz phaseChange model
         virtual void correct();
diff --git a/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.C b/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.C
index 2d576b9903..ec21192796 100644
--- a/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.C
+++ b/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.C
@@ -67,45 +67,44 @@ Foam::compressible::cavitationModels::Merkle::Merkle
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::compressible::cavitationModels::Merkle::mDotAlphal() const
+Foam::compressible::cavitationModels::Merkle::mDotcvAlphal() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
     const volScalarField::Internal mvCoeff_
     (
-        Cv_*phases_.rho1()()/(0.5*sqr(UInf_)*tInf_*phases_.rho2()())
+        Cv_*rhol()/(0.5*sqr(UInf_)*tInf_*rhov())
     );
 
     return Pair<tmp<volScalarField::Internal>>
     (
         mcCoeff_*max(p - pSat(), p0_),
-        mvCoeff_*min(p - pSat(), p0_)
+       -mvCoeff_*min(p - pSat(), p0_)
     );
 }
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::compressible::cavitationModels::Merkle::mDotP() const
+Foam::compressible::cavitationModels::Merkle::mDotcvP() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     const volScalarField::Internal mvCoeff_
     (
-        Cv_*phases_.rho1()()/(0.5*sqr(UInf_)*tInf_*phases_.rho2()())
+        Cv_*rhol()/(0.5*sqr(UInf_)*tInf_*rhov())
     );
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        mcCoeff_*(1.0 - limitedAlpha1)*pos0(p - pSat()),
-
-        (-mvCoeff_)*limitedAlpha1*neg(p - pSat())
+        mcCoeff_*(1.0 - limitedAlphal)*pos0(p - pSat()),
+        (-mvCoeff_)*limitedAlphal*neg(p - pSat())
     );
 }
 
diff --git a/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.H b/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.H
index 42538495c0..4ff2f750b2 100644
--- a/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.H
+++ b/src/twoPhaseModels/compressibleCavitation/Merkle/Merkle.H
@@ -98,13 +98,13 @@ public:
     // Member Functions
 
         //- Return the mass condensation and vaporisation rates as a
-        //  coefficient to multiply (1 - alphal) for the condensation rate
-        //  and a coefficient to multiply  alphal for the vaporisation rate
-        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const;
+        //  coefficient to multiply alphav for the condensation rate and a
+        //  coefficient to multiply alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvAlphal() const;
 
         //- Return the mass condensation and vaporisation rates as coefficients
         //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField::Internal>> mDotP() const;
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvP() const;
 
         //- Correct the Merkle phaseChange model
         virtual void correct();
diff --git a/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.C b/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.C
index 8056241eb8..8bc3965404 100644
--- a/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.C
+++ b/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.C
@@ -73,13 +73,13 @@ Foam::compressible::cavitationModels::SchnerrSauer::SchnerrSauer
 Foam::tmp<Foam::volScalarField::Internal>
 Foam::compressible::cavitationModels::SchnerrSauer::rRb
 (
-    const volScalarField::Internal& limitedAlpha1
+    const volScalarField::Internal& limitedAlphal
 ) const
 {
     return pow
     (
         ((4*constant::mathematical::pi*n_)/3)
-       *limitedAlpha1/(1.0 + alphaNuc() - limitedAlpha1),
+       *limitedAlphal/(1.0 + alphaNuc() - limitedAlphal),
         1.0/3.0
     );
 }
@@ -99,65 +99,63 @@ Foam::compressible::cavitationModels::SchnerrSauer::pCoeff
     const volScalarField::Internal& p
 ) const
 {
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     const volScalarField::Internal rho
     (
-        limitedAlpha1*phases_.rho1()()
-      + (scalar(1) - limitedAlpha1)*phases_.rho2()()
+        limitedAlphal*rhol()
+      + (scalar(1) - limitedAlphal)*rhov()
     );
 
     return
-        (3*phases_.rho1()()*phases_.rho2()())*sqrt(2/(3*phases_.rho1()()))
-       *rRb(limitedAlpha1)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));
+        (3*rhol()*rhov())*sqrt(2/(3*rhol()))
+       *rRb(limitedAlphal)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));
 }
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::compressible::cavitationModels::SchnerrSauer::mDotAlphal() const
+Foam::compressible::cavitationModels::SchnerrSauer::mDotcvAlphal() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
     const volScalarField::Internal pCoeff(this->pCoeff(p));
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        Cc_*limitedAlpha1*pCoeff*max(p - pSat(), p0_),
-
-        Cv_*(1.0 + alphaNuc() - limitedAlpha1)*pCoeff*min(p - pSat(), p0_)
+        Cc_*limitedAlphal*pCoeff*max(p - pSat(), p0_),
+       -Cv_*(1.0 + alphaNuc() - limitedAlphal)*pCoeff*min(p - pSat(), p0_)
     );
 }
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::compressible::cavitationModels::SchnerrSauer::mDotP() const
+Foam::compressible::cavitationModels::SchnerrSauer::mDotcvP() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
     const volScalarField::Internal pCoeff(this->pCoeff(p));
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
-    const volScalarField::Internal apCoeff(limitedAlpha1*pCoeff);
+    const volScalarField::Internal apCoeff(limitedAlphal*pCoeff);
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        Cc_*(1.0 - limitedAlpha1)*pos0(p - pSat())*apCoeff,
-
-        (-Cv_)*(1.0 + alphaNuc() - limitedAlpha1)*neg(p - pSat())*apCoeff
+        Cc_*(1.0 - limitedAlphal)*pos0(p - pSat())*apCoeff,
+        (-Cv_)*(1.0 + alphaNuc() - limitedAlphal)*neg(p - pSat())*apCoeff
     );
 }
 
diff --git a/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.H b/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.H
index f77c5941bd..96ff88a761 100644
--- a/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.H
+++ b/src/twoPhaseModels/compressibleCavitation/SchnerrSauer/SchnerrSauer.H
@@ -82,9 +82,9 @@ class SchnerrSauer
         dimensionedScalar alphaNuc() const;
 
         //- Reciprocal bubble radius
-        tmp<volScalarField::Internal>rRb
+        tmp<volScalarField::Internal> rRb
         (
-            const volScalarField::Internal& limitedAlpha1
+            const volScalarField::Internal& limitedAlphal
         ) const;
 
         //- Part of the condensation and vapourisation rates
@@ -118,13 +118,13 @@ public:
     // Member Functions
 
         //- Return the mass condensation and vaporisation rates as a
-        //  coefficient to multiply (1 - alphal) for the condensation rate
-        //  and a coefficient to multiply  alphal for the vaporisation rate
-        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const;
+        //  coefficient to multiply alphav for the condensation rate and a
+        //  coefficient to multiply alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvAlphal() const;
 
         //- Return the mass condensation and vaporisation rates as coefficients
         //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField::Internal>> mDotP() const;
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvP() const;
 
         //- Correct the SchnerrSauer phaseChange model
         virtual void correct();
diff --git a/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.C b/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.C
index 6b1771156d..e941a22df2 100644
--- a/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.C
+++ b/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.C
@@ -47,6 +47,7 @@ Foam::compressible::cavitationModel::cavitationModel
 )
 :
     phases_(phases),
+    liquidIndex_(phases.index(dict.lookup<word>("liquid"))),
     pSat_("pSat", dimPressure, dict.lookup("pSat"))
 {}
 
diff --git a/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.H b/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.H
index a37a84df16..19b3181b29 100644
--- a/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.H
+++ b/src/twoPhaseModels/compressibleCavitation/cavitationModel/cavitationModel.H
@@ -56,13 +56,43 @@ protected:
 
     // Protected data
 
-        //- Mixture properties
+        //- Phases
         const compressibleTwoPhases& phases_;
 
+        //- Index of the liquid
+        const bool liquidIndex_;
+
         //- Saturation vapour pressure
         dimensionedScalar pSat_;
 
 
+    // Protected Member Functions
+
+        //- Return the liquid density
+        inline const volScalarField::Internal& alphal() const
+        {
+            return phases_.alpha(liquidIndex_);
+        }
+
+        //- Return the vapour density
+        inline const volScalarField::Internal& alphav() const
+        {
+            return phases_.alpha(!liquidIndex_);
+        }
+
+        //- Return the liquid density
+        inline const volScalarField::Internal& rhol() const
+        {
+            return phases_.rho(liquidIndex_);
+        }
+
+        //- Return the vapour density
+        inline const volScalarField::Internal& rhov() const
+        {
+            return phases_.rho(!liquidIndex_);
+        }
+
+
 public:
 
     //- Runtime type information
@@ -109,20 +139,34 @@ public:
 
     // Member Functions
 
-        //- Return const-access to the saturation vapour pressure
-        const dimensionedScalar& pSat() const
+        //- Return the saturation vapour pressure
+        inline const dimensionedScalar& pSat() const
         {
             return pSat_;
         }
 
         //- Return the mass condensation and vaporisation rates as a
-        //  coefficient to multiply (1 - alphal) for the condensation rate
-        //  and a coefficient to multiply  alphal for the vaporisation rate
-        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const = 0;
+        //  coefficient to multiply alphav for the condensation rate and a
+        //  coefficient to multiply alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvAlphal() const = 0;
 
         //- Return the mass condensation and vaporisation rates as coefficients
         //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField::Internal>> mDotP() const = 0;
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvP() const = 0;
+
+        //- Return the mass transfer rates of the two phases as coefficients to
+        //  multiply the volume fraction of the other phase
+        inline Pair<tmp<volScalarField::Internal>> mDot12Alpha() const
+        {
+            return liquidIndex_ ? reverse(mDotcvAlphal()) : mDotcvAlphal();
+        }
+
+        //- Return the mass transfer rates of the two phases as coefficients to
+        //  multiply (p - pSat)
+        inline Pair<tmp<volScalarField::Internal>> mDot12P() const
+        {
+            return liquidIndex_ ? reverse(mDotcvP()) : mDotcvP();
+        }
 
         //- Correct the cavitation model
         virtual void correct() = 0;
diff --git a/src/twoPhaseModels/compressibleTwoPhaseMixture/compressibleTwoPhases.H b/src/twoPhaseModels/compressibleTwoPhaseMixture/compressibleTwoPhases.H
index afc3072c52..d908b9a45a 100644
--- a/src/twoPhaseModels/compressibleTwoPhaseMixture/compressibleTwoPhases.H
+++ b/src/twoPhaseModels/compressibleTwoPhaseMixture/compressibleTwoPhases.H
@@ -77,11 +77,23 @@ public:
         //- Return the thermo for phase 2
         virtual const rhoThermo& thermo2() const = 0;
 
+        //- Return the density of a given phase
+        inline const rhoThermo& thermo(const bool index) const
+        {
+            return index ? thermo2() : thermo1();
+        }
+
         //- Return the density of phase 1
         virtual const volScalarField& rho1() const = 0;
 
         //- Return the density of phase 2
         virtual const volScalarField& rho2() const = 0;
+
+        //- Return the density of a given phase
+        inline const volScalarField& rho(const bool index) const
+        {
+            return index ? rho2() : rho1();
+        }
 };
 
 
diff --git a/src/twoPhaseModels/incompressibleCavitation/Kunz/Kunz.C b/src/twoPhaseModels/incompressibleCavitation/Kunz/Kunz.C
index 0e55c6f6be..4ee67be2bb 100644
--- a/src/twoPhaseModels/incompressibleCavitation/Kunz/Kunz.C
+++ b/src/twoPhaseModels/incompressibleCavitation/Kunz/Kunz.C
@@ -54,8 +54,8 @@ Foam::cavitationModels::Kunz::Kunz
 
     p0_("0", pSat().dimensions(), 0.0),
 
-    mcCoeff_(Cc_*phases_.rho2()/tInf_),
-    mvCoeff_(Cv_*phases_.rho2()/(0.5*phases_.rho1()*sqr(UInf_)*tInf_))
+    mcCoeff_(Cc_*rhov()/tInf_),
+    mvCoeff_(Cv_*rhov()/(0.5*rhol()*sqr(UInf_)*tInf_))
 {
     correct();
 }
@@ -64,43 +64,41 @@ Foam::cavitationModels::Kunz::Kunz
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::cavitationModels::Kunz::mDotAlphal() const
+Foam::cavitationModels::Kunz::mDotcvAlpha() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        mcCoeff_*sqr(limitedAlpha1)
+        mcCoeff_*sqr(limitedAlphal)
        *max(p - pSat(), p0_)/max(p - pSat(), 0.01*pSat()),
-
-        mvCoeff_*min(p - pSat(), p0_)
+       -mvCoeff_*min(p - pSat(), p0_)
     );
 }
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::cavitationModels::Kunz::mDotP() const
+Foam::cavitationModels::Kunz::mDotcvP() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        mcCoeff_*sqr(limitedAlpha1)*(1.0 - limitedAlpha1)
+        mcCoeff_*sqr(limitedAlphal)*(1.0 - limitedAlphal)
        *pos0(p - pSat())/max(p - pSat(), 0.01*pSat()),
-
-        (-mvCoeff_)*limitedAlpha1*neg(p - pSat())
+        (-mvCoeff_)*limitedAlphal*neg(p - pSat())
     );
 }
 
@@ -118,8 +116,8 @@ bool Foam::cavitationModels::Kunz::read(const dictionary& dict)
         dict.lookup("Cc") >> Cc_;
         dict.lookup("Cv") >> Cv_;
 
-        mcCoeff_ = Cc_*phases_.rho2()/tInf_;
-        mvCoeff_ = Cv_*phases_.rho2()/(0.5*phases_.rho1()*sqr(UInf_)*tInf_);
+        mcCoeff_ = Cc_*rhov()/tInf_;
+        mvCoeff_ = Cv_*rhov()/(0.5*rhol()*sqr(UInf_)*tInf_);
 
         return true;
     }
diff --git a/src/twoPhaseModels/incompressibleCavitation/Kunz/Kunz.H b/src/twoPhaseModels/incompressibleCavitation/Kunz/Kunz.H
index 1853bdd066..d346c9c1c7 100644
--- a/src/twoPhaseModels/incompressibleCavitation/Kunz/Kunz.H
+++ b/src/twoPhaseModels/incompressibleCavitation/Kunz/Kunz.H
@@ -103,13 +103,13 @@ public:
     // Member Functions
 
         //- Return the mass condensation and vaporisation rates as a
-        //  coefficient to multiply (1 - alphal) for the condensation rate
-        //  and a coefficient to multiply  alphal for the vaporisation rate
-        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const;
+        //  coefficient to multiply alphav for the condensation rate and a
+        //  coefficient to multiply alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvAlpha() const;
 
         //- Return the mass condensation and vaporisation rates as coefficients
         //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField::Internal>> mDotP() const;
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvP() const;
 
         //- Correct the Kunz phaseChange model
         virtual void correct();
diff --git a/src/twoPhaseModels/incompressibleCavitation/Merkle/Merkle.C b/src/twoPhaseModels/incompressibleCavitation/Merkle/Merkle.C
index d8a1650db8..493fc027e2 100644
--- a/src/twoPhaseModels/incompressibleCavitation/Merkle/Merkle.C
+++ b/src/twoPhaseModels/incompressibleCavitation/Merkle/Merkle.C
@@ -55,7 +55,7 @@ Foam::cavitationModels::Merkle::Merkle
     p0_("0", pSat().dimensions(), 0.0),
 
     mcCoeff_(Cc_/(0.5*sqr(UInf_)*tInf_)),
-    mvCoeff_(Cv_*phases_.rho1()/(0.5*sqr(UInf_)*tInf_*phases_.rho2()))
+    mvCoeff_(Cv_*rhol()/(0.5*sqr(UInf_)*tInf_*rhov()))
 {
     correct();
 }
@@ -64,7 +64,7 @@ Foam::cavitationModels::Merkle::Merkle
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::cavitationModels::Merkle::mDotAlphal() const
+Foam::cavitationModels::Merkle::mDotcvAlpha() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
@@ -72,26 +72,26 @@ Foam::cavitationModels::Merkle::mDotAlphal() const
     return Pair<tmp<volScalarField::Internal>>
     (
         mcCoeff_*max(p - pSat(), p0_),
-        mvCoeff_*min(p - pSat(), p0_)
+       -mvCoeff_*min(p - pSat(), p0_)
     );
 }
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::cavitationModels::Merkle::mDotP() const
+Foam::cavitationModels::Merkle::mDotcvP() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        mcCoeff_*(1.0 - limitedAlpha1)*pos0(p - pSat()),
-        (-mvCoeff_)*limitedAlpha1*neg(p - pSat())
+        mcCoeff_*(1.0 - limitedAlphal)*pos0(p - pSat()),
+        (-mvCoeff_)*limitedAlphal*neg(p - pSat())
     );
 }
 
@@ -110,7 +110,7 @@ bool Foam::cavitationModels::Merkle::read(const dictionary& dict)
         dict.lookup("Cv") >> Cv_;
 
         mcCoeff_ = Cc_/(0.5*sqr(UInf_)*tInf_);
-        mvCoeff_ = Cv_*phases_.rho1()/(0.5*sqr(UInf_)*tInf_*phases_.rho2());
+        mvCoeff_ = Cv_*rhol()/(0.5*sqr(UInf_)*tInf_*rhov());
 
         return true;
     }
diff --git a/src/twoPhaseModels/incompressibleCavitation/Merkle/Merkle.H b/src/twoPhaseModels/incompressibleCavitation/Merkle/Merkle.H
index 5a815de952..bc348afbf3 100644
--- a/src/twoPhaseModels/incompressibleCavitation/Merkle/Merkle.H
+++ b/src/twoPhaseModels/incompressibleCavitation/Merkle/Merkle.H
@@ -97,13 +97,13 @@ public:
     // Member Functions
 
         //- Return the mass condensation and vaporisation rates as a
-        //  coefficient to multiply (1 - alphal) for the condensation rate
-        //  and a coefficient to multiply  alphal for the vaporisation rate
-        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const;
+        //  coefficient to multiply alphav for the condensation rate and a
+        //  coefficient to multiply alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvAlpha() const;
 
         //- Return the mass condensation and vaporisation rates as coefficients
         //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField::Internal>> mDotP() const;
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvP() const;
 
         //- Correct the Merkle phaseChange model
         virtual void correct();
diff --git a/src/twoPhaseModels/incompressibleCavitation/SchnerrSauer/SchnerrSauer.C b/src/twoPhaseModels/incompressibleCavitation/SchnerrSauer/SchnerrSauer.C
index 76dd22c42e..dc2cc7ae86 100644
--- a/src/twoPhaseModels/incompressibleCavitation/SchnerrSauer/SchnerrSauer.C
+++ b/src/twoPhaseModels/incompressibleCavitation/SchnerrSauer/SchnerrSauer.C
@@ -70,13 +70,13 @@ Foam::cavitationModels::SchnerrSauer::SchnerrSauer
 Foam::tmp<Foam::volScalarField::Internal>
 Foam::cavitationModels::SchnerrSauer::rRb
 (
-    const volScalarField::Internal& limitedAlpha1
+    const volScalarField::Internal& limitedAlphal
 ) const
 {
     return pow
     (
         ((4*constant::mathematical::pi*n_)/3)
-       *limitedAlpha1/(1.0 + alphaNuc() - limitedAlpha1),
+       *limitedAlphal/(1.0 + alphaNuc() - limitedAlphal),
         1.0/3.0
     );
 }
@@ -96,65 +96,63 @@ Foam::cavitationModels::SchnerrSauer::pCoeff
     const volScalarField::Internal& p
 ) const
 {
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     const volScalarField::Internal rho
     (
-        limitedAlpha1*phases_.rho1()
-      + (scalar(1) - limitedAlpha1)*phases_.rho2()
+        limitedAlphal*rhol()
+      + (scalar(1) - limitedAlphal)*rhov()
     );
 
     return
-        (3*phases_.rho1()*phases_.rho2())*sqrt(2/(3*phases_.rho1()))
-       *rRb(limitedAlpha1)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));
+        (3*rhol()*rhov())*sqrt(2/(3*rhol()))
+       *rRb(limitedAlphal)/(rho*sqrt(mag(p - pSat()) + 0.01*pSat()));
 }
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::cavitationModels::SchnerrSauer::mDotAlphal() const
+Foam::cavitationModels::SchnerrSauer::mDotcvAlpha() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
     const volScalarField::Internal pCoeff(this->pCoeff(p));
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        Cc_*limitedAlpha1*pCoeff*max(p - pSat(), p0_),
-
-        Cv_*(1.0 + alphaNuc() - limitedAlpha1)*pCoeff*min(p - pSat(), p0_)
+        Cc_*limitedAlphal*pCoeff*max(p - pSat(), p0_),
+       -Cv_*(1.0 + alphaNuc() - limitedAlphal)*pCoeff*min(p - pSat(), p0_)
     );
 }
 
 
 Foam::Pair<Foam::tmp<Foam::volScalarField::Internal>>
-Foam::cavitationModels::SchnerrSauer::mDotP() const
+Foam::cavitationModels::SchnerrSauer::mDotcvP() const
 {
     const volScalarField::Internal& p =
         phases_.mesh().lookupObject<volScalarField>("p");
 
     const volScalarField::Internal pCoeff(this->pCoeff(p));
 
-    const volScalarField::Internal limitedAlpha1
+    const volScalarField::Internal limitedAlphal
     (
-        min(max(phases_.alpha1()(), scalar(0)), scalar(1))
+        min(max(alphal(), scalar(0)), scalar(1))
     );
 
-    const volScalarField::Internal apCoeff(limitedAlpha1*pCoeff);
+    const volScalarField::Internal apCoeff(limitedAlphal*pCoeff);
 
     return Pair<tmp<volScalarField::Internal>>
     (
-        Cc_*(1.0 - limitedAlpha1)*pos0(p - pSat())*apCoeff,
-
-        (-Cv_)*(1.0 + alphaNuc() - limitedAlpha1)*neg(p - pSat())*apCoeff
+        Cc_*(1.0 - limitedAlphal)*pos0(p - pSat())*apCoeff,
+        (-Cv_)*(1.0 + alphaNuc() - limitedAlphal)*neg(p - pSat())*apCoeff
     );
 }
 
diff --git a/src/twoPhaseModels/incompressibleCavitation/SchnerrSauer/SchnerrSauer.H b/src/twoPhaseModels/incompressibleCavitation/SchnerrSauer/SchnerrSauer.H
index 3cb70e2010..c670edaeaa 100644
--- a/src/twoPhaseModels/incompressibleCavitation/SchnerrSauer/SchnerrSauer.H
+++ b/src/twoPhaseModels/incompressibleCavitation/SchnerrSauer/SchnerrSauer.H
@@ -80,9 +80,9 @@ class SchnerrSauer
         dimensionedScalar alphaNuc() const;
 
         //- Reciprocal bubble radius
-        tmp<volScalarField::Internal>rRb
+        tmp<volScalarField::Internal> rRb
         (
-            const volScalarField::Internal& limitedAlpha1
+            const volScalarField::Internal& limitedAlphal
         ) const;
 
         //- Part of the condensation and vapourisation rates
@@ -116,13 +116,13 @@ public:
     // Member Functions
 
         //- Return the mass condensation and vaporisation rates as a
-        //  coefficient to multiply (1 - alphal) for the condensation rate
-        //  and a coefficient to multiply  alphal for the vaporisation rate
-        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const;
+        //  coefficient to multiply alphav for the condensation rate and a
+        //  coefficient to multiply alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvAlpha() const;
 
         //- Return the mass condensation and vaporisation rates as coefficients
         //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField::Internal>> mDotP() const;
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvP() const;
 
         //- Correct the SchnerrSauer phaseChange model
         virtual void correct();
diff --git a/src/twoPhaseModels/incompressibleCavitation/cavitationModel/cavitationModel.C b/src/twoPhaseModels/incompressibleCavitation/cavitationModel/cavitationModel.C
index d2804e7636..b6767f88be 100644
--- a/src/twoPhaseModels/incompressibleCavitation/cavitationModel/cavitationModel.C
+++ b/src/twoPhaseModels/incompressibleCavitation/cavitationModel/cavitationModel.C
@@ -44,6 +44,7 @@ Foam::cavitationModel::cavitationModel
 )
 :
     phases_(phases),
+    liquidIndex_(phases.index(dict.lookup<word>("liquid"))),
     pSat_("pSat", dimPressure, dict.lookup("pSat"))
 {}
 
diff --git a/src/twoPhaseModels/incompressibleCavitation/cavitationModel/cavitationModel.H b/src/twoPhaseModels/incompressibleCavitation/cavitationModel/cavitationModel.H
index 1e4bd821df..91981fb2cc 100644
--- a/src/twoPhaseModels/incompressibleCavitation/cavitationModel/cavitationModel.H
+++ b/src/twoPhaseModels/incompressibleCavitation/cavitationModel/cavitationModel.H
@@ -57,10 +57,40 @@ protected:
         //- Mixture properties
         const incompressibleTwoPhases& phases_;
 
+        //- Index of the liquid
+        const bool liquidIndex_;
+
         //- Saturation vapour pressure
         dimensionedScalar pSat_;
 
 
+    // Protected Member Functions
+
+        //- Return the liquid density
+        inline const volScalarField::Internal& alphal() const
+        {
+            return phases_.alpha(liquidIndex_);
+        }
+
+        //- Return the vapour density
+        inline const volScalarField::Internal& alphav() const
+        {
+            return phases_.alpha(!liquidIndex_);
+        }
+
+        //- Return the liquid density
+        inline const dimensionedScalar& rhol() const
+        {
+            return phases_.rho(liquidIndex_);
+        }
+
+        //- Return the vapour density
+        inline const dimensionedScalar& rhov() const
+        {
+            return phases_.rho(!liquidIndex_);
+        }
+
+
 public:
 
     //- Runtime type information
@@ -107,20 +137,34 @@ public:
 
     // Member Functions
 
-        //- Return const-access to the saturation vapour pressure
-        const dimensionedScalar& pSat() const
+        //- Return the saturation vapour pressure
+        inline const dimensionedScalar& pSat() const
         {
             return pSat_;
         }
 
         //- Return the mass condensation and vaporisation rates as a
-        //  coefficient to multiply (1 - alphal) for the condensation rate
-        //  and a coefficient to multiply  alphal for the vaporisation rate
-        virtual Pair<tmp<volScalarField::Internal>> mDotAlphal() const = 0;
+        //  coefficient to multiply alphav for the condensation rate and a
+        //  coefficient to multiply alphal for the vaporisation rate
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvAlpha() const = 0;
 
         //- Return the mass condensation and vaporisation rates as coefficients
         //  to multiply (p - pSat)
-        virtual Pair<tmp<volScalarField::Internal>> mDotP() const = 0;
+        virtual Pair<tmp<volScalarField::Internal>> mDotcvP() const = 0;
+
+        //- Return the mass transfer rates of the two phases as coefficients to
+        //  multiply the volume fraction of the other phase
+        inline Pair<tmp<volScalarField::Internal>> mDot12Alpha() const
+        {
+            return liquidIndex_ ? reverse(mDotcvAlpha()) : mDotcvAlpha();
+        }
+
+        //- Return the mass transfer rates of the two phases as coefficients to
+        //  multiply (p - pSat)
+        inline Pair<tmp<volScalarField::Internal>> mDot12P() const
+        {
+            return liquidIndex_ ? reverse(mDotcvP()) : mDotcvP();
+        }
 
         //- Correct the cavitation model
         virtual void correct() = 0;
diff --git a/src/twoPhaseModels/incompressibleTwoPhaseMixture/incompressibleTwoPhases.H b/src/twoPhaseModels/incompressibleTwoPhaseMixture/incompressibleTwoPhases.H
index af92b75842..87d1b6464d 100644
--- a/src/twoPhaseModels/incompressibleTwoPhaseMixture/incompressibleTwoPhases.H
+++ b/src/twoPhaseModels/incompressibleTwoPhaseMixture/incompressibleTwoPhases.H
@@ -52,6 +52,7 @@ class incompressibleTwoPhases
 {
 public:
 
+    //- Runtime type information
     TypeName("incompressibleTwoPhases");
 
 
@@ -69,11 +70,17 @@ public:
 
     // Member Functions
 
-        //- Return const-access to phase1 density
+        //- Return the density of phase1
         virtual const dimensionedScalar& rho1() const = 0;
 
-        //- Return const-access to phase2 density
+        //- Return the density of phase2
         virtual const dimensionedScalar& rho2() const = 0;
+
+        //- Return the density of a given phase
+        inline const dimensionedScalar& rho(const bool index) const
+        {
+            return index ? rho2() : rho1();
+        }
 };
 
 
diff --git a/src/twoPhaseModels/twoPhaseMixture/twoPhases.H b/src/twoPhaseModels/twoPhaseMixture/twoPhases.H
index bbd71ecd42..a1b16add33 100644
--- a/src/twoPhaseModels/twoPhaseMixture/twoPhases.H
+++ b/src/twoPhaseModels/twoPhaseMixture/twoPhases.H
@@ -50,6 +50,7 @@ class twoPhases
 {
 public:
 
+    //- Runtime type information
     TypeName("twoPhases");
 
 
@@ -73,11 +74,37 @@ public:
             return alpha1().mesh();
         }
 
-        //- Return the phase-fraction of phase 1
+        //- Return the index associated with a given phase name
+        inline bool index(const word& name) const
+        {
+            if (alpha1().group() == name)
+            {
+                return false;
+            }
+
+            if (alpha2().group() == name)
+            {
+                return true;
+            }
+
+            FatalErrorInFunction
+                << "Index of phase " << name << " requested from phases"
+                << alpha1().group() << " and " << alpha2().group()
+                << exit(FatalError);
+            return false;
+        }
+
+        //- Return the volume-fraction of phase 1
         virtual const volScalarField& alpha1() const = 0;
 
-        //- Return the phase-fraction of phase 2
+        //- Return the volume-fraction of phase 2
         virtual const volScalarField& alpha2() const = 0;
+
+        //- Return the volume-fraction of a given phase
+        const volScalarField& alpha(const bool index) const
+        {
+            return index ? alpha2() : alpha1();
+        }
 };
 
 
diff --git a/tutorials/modules/compressibleVoF/ballValve/0/alpha.water b/tutorials/modules/compressibleVoF/ballValve/0/alpha.vapour
similarity index 95%
rename from tutorials/modules/compressibleVoF/ballValve/0/alpha.water
rename to tutorials/modules/compressibleVoF/ballValve/0/alpha.vapour
index b544f8a1e7..781ff13fc3 100644
--- a/tutorials/modules/compressibleVoF/ballValve/0/alpha.water
+++ b/tutorials/modules/compressibleVoF/ballValve/0/alpha.vapour
@@ -10,13 +10,13 @@ FoamFile
     format      ascii;
     class       volScalarField;
     location    "0";
-    object      alpha.water;
+    object      alpha.vapour;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [0 0 0 0 0 0 0];
 
-internalField   uniform 1;
+internalField   uniform 0;
 
 boundaryField
 {
diff --git a/tutorials/modules/compressibleVoF/ballValve/constant/fvModels b/tutorials/modules/compressibleVoF/ballValve/constant/fvModels
index e4daa01fa8..5606fa6bfe 100644
--- a/tutorials/modules/compressibleVoF/ballValve/constant/fvModels
+++ b/tutorials/modules/compressibleVoF/ballValve/constant/fvModels
@@ -22,6 +22,8 @@ VoFCavitation
 
     model   SchnerrSauer;
 
+    liquid  water;
+
     pSat    2300;
 
     n       1.6e+13;
diff --git a/tutorials/modules/compressibleVoF/ballValve/constant/phaseProperties b/tutorials/modules/compressibleVoF/ballValve/constant/phaseProperties
index 1fbc9cd412..cc486a44ef 100644
--- a/tutorials/modules/compressibleVoF/ballValve/constant/phaseProperties
+++ b/tutorials/modules/compressibleVoF/ballValve/constant/phaseProperties
@@ -14,7 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-phases          (water vapour);
+phases          (vapour water);
 
 pMin            100;
 
diff --git a/tutorials/modules/compressibleVoF/ballValve/system/fvSolution b/tutorials/modules/compressibleVoF/ballValve/system/fvSolution
index 525a9df5a8..64ff01f630 100644
--- a/tutorials/modules/compressibleVoF/ballValve/system/fvSolution
+++ b/tutorials/modules/compressibleVoF/ballValve/system/fvSolution
@@ -16,7 +16,7 @@ FoamFile
 
 solvers
 {
-    "alpha.water.*"
+    "alpha.vapour.*"
     {
         nAlphaCorr      2;
         nAlphaSubCycles 1;
diff --git a/tutorials/multiphase/interFoam/RAS/propeller/constant/fvModels b/tutorials/multiphase/interFoam/RAS/propeller/constant/fvModels
index 922d18b16f..d2a9b6bb02 100644
--- a/tutorials/multiphase/interFoam/RAS/propeller/constant/fvModels
+++ b/tutorials/multiphase/interFoam/RAS/propeller/constant/fvModels
@@ -26,6 +26,8 @@ VoFCavitation
     {
         pSat    2300;   // Saturation pressure
 
+        liquid  water;
+
         UInf    20.0;
         tInf    0.005; // L = 0.1 m
         Cc      1000;
@@ -36,6 +38,8 @@ VoFCavitation
     {
         pSat    2300;   // Saturation pressure
 
+        liquid  water;
+
         UInf    20.0;
         tInf    0.005;  // L = 0.1 m
         Cc      80;
@@ -46,6 +50,8 @@ VoFCavitation
     {
         pSat    2300;   // Saturation pressure
 
+        liquid  water;
+
         n       1.6e+13;
         dNuc    2.0e-06;
         Cc      1;
diff --git a/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/fvModels b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/fvModels
index 922d18b16f..d2a9b6bb02 100644
--- a/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/fvModels
+++ b/tutorials/multiphase/interFoam/laminar/cavitatingBullet/constant/fvModels
@@ -26,6 +26,8 @@ VoFCavitation
     {
         pSat    2300;   // Saturation pressure
 
+        liquid  water;
+
         UInf    20.0;
         tInf    0.005; // L = 0.1 m
         Cc      1000;
@@ -36,6 +38,8 @@ VoFCavitation
     {
         pSat    2300;   // Saturation pressure
 
+        liquid  water;
+
         UInf    20.0;
         tInf    0.005;  // L = 0.1 m
         Cc      80;
@@ -46,6 +50,8 @@ VoFCavitation
     {
         pSat    2300;   // Saturation pressure
 
+        liquid  water;
+
         n       1.6e+13;
         dNuc    2.0e-06;
         Cc      1;