diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C
index 4b52881506..28ac209b9c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/diameterModels/velocityGroup/sizeGroup/shapeModels/fractal/fractal.C
@@ -210,8 +210,16 @@ void Foam::diameterModels::shapeModels::fractal::correct()
       - sinteringModel_->R()
       + Su_
       - fvm::Sp(popBal.Sp(fi.i()())*fi, kappa_)
-      + fvc::ddt(fi.phase().residualAlpha(), kappa_)
-      - fvm::ddt(fi.phase().residualAlpha(), kappa_)
+
+      - correction
+        (
+            fvm::Sp
+            (
+                max(phase.residualAlpha() - alpha, scalar(0))
+               /sizeGroup_.mesh().time().deltaT(),
+                kappa_
+            )
+        )
     );
 
     kappaEqn.relax();
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
index a614b53918..189a45a4fa 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/MultiComponentPhaseModel/MultiComponentPhaseModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -110,8 +110,15 @@ Foam::MultiComponentPhaseModel<BasePhaseModel>::YiEqn(volScalarField& Yi)
      ==
         alpha*this->R(Yi)
 
-      + fvc::ddt(residualAlpha_*rho, Yi)
-      - fvm::ddt(residualAlpha_*rho, Yi)
+      - correction
+        (
+            fvm::Sp
+            (
+                max(residualAlpha_ - alpha, scalar(0))*rho
+               /this->mesh().time().deltaT(),
+                Yi
+            )
+        )
     );
 }
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
index adc6f8eb46..5261cca161 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
@@ -1066,7 +1066,6 @@ void Foam::diameterModels::populationBalanceModel::solve()
                 sizeGroup& fi = sizeGroups_[i];
                 const phaseModel& phase = fi.phase();
                 const volScalarField& alpha = phase;
-                const dimensionedScalar& residualAlpha = phase.residualAlpha();
                 const volScalarField& rho = phase.thermo().rho();
 
                 fvScalarMatrix sizeGroupEqn
@@ -1077,8 +1076,16 @@ void Foam::diameterModels::populationBalanceModel::solve()
                     Su_[i]
                   - fvm::Sp(Sp_[i], fi)
                   + fluid_.fvModels().source(alpha, rho, fi)/rho
-                  + fvc::ddt(residualAlpha, fi)
-                  - fvm::ddt(residualAlpha, fi)
+
+                  - correction
+                    (
+                        fvm::Sp
+                        (
+                            max(phase.residualAlpha() - alpha, scalar(0))
+                           /this->mesh().time().deltaT(),
+                            fi
+                        )
+                    )
                 );
 
                 sizeGroupEqn.relax();