diff --git a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
index b7608d1d0d..c7fbe09556 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H
@@ -1,34 +1,33 @@
-    fvVectorMatrix UEqn
+fvVectorMatrix UEqn
+(
+    fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - fvm::Sp(contErr, U)
+  + MRF.DDt(rho, U)
+  + turbulence.divDevTau(U)
+ ==
+    fvModels.source(rho, U)
+);
+
+UEqn.relax();
+
+fvConstraints.constrain(UEqn);
+
+if (pimple.momentumPredictor())
+{
+    solve
     (
-        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      - fvm::Sp(contErr, U)
-      + MRF.DDt(rho, U)
-      + turbulence.divDevTau(U)
+        UEqn
      ==
-        fvModels.source(rho, U)
+        fvc::reconstruct
+        (
+            (
+                mixture.surfaceTensionForce()
+              - ghf*fvc::snGrad(rho)
+              - fvc::snGrad(p_rgh)
+            ) * mesh.magSf()
+        )
     );
 
-    UEqn.relax();
+    fvConstraints.constrain(U);
 
-    fvConstraints.constrain(UEqn);
-
-    if (pimple.momentumPredictor())
-    {
-        solve
-        (
-            UEqn
-         ==
-            fvc::reconstruct
-            (
-                (
-                    mixture.surfaceTensionForce()
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                ) * mesh.magSf()
-            )
-        );
-
-        fvConstraints.constrain(U);
-
-        K = 0.5*magSqr(U);
-    }
+    K = 0.5*magSqr(U);
+}
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H
index e4fef9717e..27423fb975 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleAlphaEqnSubCycle.H
@@ -61,7 +61,7 @@ volScalarField::Internal contErr
 (
     (
         fvc::ddt(rho) + fvc::div(rhoPhi)
-      - (fvModels.source(alpha1, mixture.thermo1().rho())&rho1)
-      - (fvModels.source(alpha2, mixture.thermo2().rho())&rho2)
+      - (fvModels.source(alpha1, rho1)&rho1)
+      - (fvModels.source(alpha2, rho2)&rho2)
     )()
 );
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
index d3e8b59a86..7574df4830 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
@@ -19,4 +19,6 @@ EXE_LIBS = \
     -lmomentumTransportModels \
     -lcompressibleMomentumTransportModels \
     -lfiniteVolume \
+    -lfvModels \
+    -lfvConstraints \
     -lmeshTools
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H
index 89c2c0996a..0c42e7404f 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H
@@ -1,18 +1,25 @@
 {
     fvScalarMatrix TEqn
     (
-        fvm::ddt(rho, T)
-      + fvm::div(mixture.rhoPhi(), T)
+        fvm::ddt(rho, T) + fvm::div(mixture.rhoPhi(), T) - fvm::Sp(contErr, T)
       - fvm::laplacian(mixture.alphaEff(rho*turbulence->nut()), T)
       + (
-            fvc::div(fvc::absolute(phi, U), p)
-          + fvc::ddt(rho, K) + fvc::div(mixture.rhoPhi(), K)
-        )*mixture.rCv()
+            fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p
+          + (fvc::ddt(rho, K) + fvc::div(mixture.rhoPhi(), K))()()
+          - (U()&(fvModels.source(rho, U)&U)()) - contErr*K
+        )*mixture.rCv()()
+     ==
+        fvModels.source(rho, T)
     );
 
     TEqn.relax();
+
+    fvConstraints.constrain(TEqn);
+
     TEqn.solve();
 
+    fvConstraints.constrain(T);
+
     mixture.correctThermo();
     mixture.correct();
 }
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H
index ed987ced7a..936820c8a4 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H
@@ -1,27 +1,32 @@
-    fvVectorMatrix UEqn
+fvVectorMatrix UEqn
+(
+    fvm::ddt(rho, U) + fvm::div(mixture.rhoPhi(), U) - fvm::Sp(contErr, U)
+  + turbulence->divDevTau(U)
+ ==
+    fvModels.source(rho, U)
+);
+
+UEqn.relax();
+
+fvConstraints.constrain(UEqn);
+
+if (pimple.momentumPredictor())
+{
+    solve
     (
-        fvm::ddt(rho, U)
-      + fvm::div(mixture.rhoPhi(), U)
-      + turbulence->divDevTau(U)
+        UEqn
+     ==
+        fvc::reconstruct
+        (
+            (
+                mixture.surfaceTensionForce()
+              - ghf*fvc::snGrad(rho)
+              - fvc::snGrad(p_rgh)
+            ) * mesh.magSf()
+        )
     );
 
-    UEqn.relax();
+    fvConstraints.constrain(U);
 
-    if (pimple.momentumPredictor())
-    {
-        solve
-        (
-            UEqn
-         ==
-            fvc::reconstruct
-            (
-                (
-                    mixture.surfaceTensionForce()
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(p_rgh)
-                ) * mesh.magSf()
-            )
-        );
-
-        K = 0.5*magSqr(U);
-    }
+    K = 0.5*magSqr(U);
+}
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
index 69429294c9..6158acb6f5 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,6 +39,8 @@ Description
 #include "compressibleMultiphaseMixture.H"
 #include "compressibleMomentumTransportModels.H"
 #include "pimpleControl.H"
+#include "fvModels.H"
+#include "fvConstraints.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,7 +82,7 @@ int main(int argc, char *argv[])
         {
             mixture.solve();
 
-            solve(fvm::ddt(rho) + fvc::div(mixture.rhoPhi()));
+            #include "contErr.H"
 
             #include "UEqn.H"
             #include "TEqn.H"
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
index ccf94010f9..c8b8b194ef 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H
@@ -57,3 +57,6 @@ autoPtr<compressible::momentumTransportModel> turbulence
 
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
+
+#include "createFvModels.H"
+#include "createFvConstraints.H"
diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
index fc3c1d017b..14f45fa30d 100644
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H
@@ -41,6 +41,7 @@
             (
                 fvc::ddt(rho) + thermo.psi()*correction(fvm::ddt(p_rgh))
               + fvc::div(phi, rho) - fvc::Sp(fvc::div(phi), rho)
+              - (fvModels.source(phase(), rho)&rho)
             ).ptr()
         );
 
@@ -113,6 +114,7 @@
             U = HbyA
               + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
             U.correctBoundaryConditions();
+            fvConstraints.constrain(U);
         }
     }