diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
index d2c02f836e..171714ef6c 100644
--- a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
+++ b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
@@ -168,37 +168,32 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                 const volScalarField& alpha = phase;
                 const label phasei = phase.index();
 
-                HbyAs.set
+                const volVectorField H
                 (
-                    phasei,
-                    constrainHbyA
+                    constrainH
                     (
-                        rAUs[phasei]
-                       *(
-                            UEqns[phasei].H()
-                          + byDt
-                            (
-                                max(phase.residualAlpha() - alpha, scalar(0))
-                               *phase.rho()
-                            )
-                           *phase.U()().oldTime()
-                        ),
+                        UEqns[phasei].H()
+                      + byDt
+                        (
+                            max(phase.residualAlpha() - alpha, scalar(0))
+                           *phase.rho()
+                        )
+                       *phase.U()().oldTime(),
+                        rAUs[phasei],
                         phase.U(),
                         p_rgh
                     )
                 );
 
+                HbyAs.set(phasei, rAUs[phasei]*H);
+
                 phiHbyAs.set
                 (
                     phasei,
                     new surfaceScalarField
                     (
                         IOobject::groupName("phiHbyA", phase.name()),
-                        rAUfs[phasei]
-                       *(
-                           fvc::flux(HbyAs[phasei]/rAUs[phasei])
-                         + ddtCorrs[phasei]
-                        )
+                        rAUfs[phasei]*(fvc::flux(H) + ddtCorrs[phasei])
                       - alpharAUfs[phasei]*phigFs[phasei]
                       - rAUfs[phasei]*Fs[phasei]
                     )
diff --git a/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C b/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C
index 3506704352..ce59f4a45c 100644
--- a/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C
+++ b/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,6 +75,52 @@ Foam::tmp<Foam::volVectorField> Foam::constrainHbyA
 }
 
 
+Foam::tmp<Foam::volVectorField> Foam::constrainH
+(
+    const tmp<volVectorField>& tH,
+    const volScalarField& rA,
+    const volVectorField& U,
+    const volScalarField& p
+)
+{
+    tmp<volVectorField> tHNew;
+
+    if (tH.isTmp())
+    {
+        tHNew = tH;
+        tHNew.ref().rename(IOobject::groupName("H", U.group()));
+    }
+    else
+    {
+        tHNew = volVectorField::New
+        (
+            IOobject::groupName("H", U.group()),
+            tH
+        );
+    }
+
+    volVectorField& H = tHNew.ref();
+    volVectorField::Boundary& Hbf = H.boundaryFieldRef();
+
+    forAll(U.boundaryField(), patchi)
+    {
+        if
+        (
+           !U.boundaryField()[patchi].assignable()
+        && !isA<fixedFluxExtrapolatedPressureFvPatchScalarField>
+            (
+                p.boundaryField()[patchi]
+            )
+        )
+        {
+            Hbf[patchi] = U.boundaryField()[patchi]/rA.boundaryField()[patchi];
+        }
+    }
+
+    return tHNew;
+}
+
+
 Foam::tmp<Foam::surfaceScalarField> Foam::constrainPhiHbyA
 (
     const tmp<surfaceScalarField>& tphiHbyA,
diff --git a/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.H b/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.H
index 90d6b0469c..f569bbfb74 100644
--- a/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.H
+++ b/src/finiteVolume/cfdTools/general/constrainHbyA/constrainHbyA.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,6 +51,14 @@ tmp<volVectorField> constrainHbyA
     const volScalarField& p
 );
 
+tmp<volVectorField> constrainH
+(
+    const tmp<volVectorField>& tH,
+    const volScalarField& rA,
+    const volVectorField& U,
+    const volScalarField& p
+);
+
 tmp<surfaceScalarField> constrainPhiHbyA
 (
     const tmp<surfaceScalarField>& tphiHbyA,