diff --git a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
index 984b52f6..8f1c2f7c 100644
--- a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
+++ b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
@@ -133,7 +133,6 @@ int main(int argc, char *argv[])
             {
                 // besides this pEqn, OF offers a "pimple consistent"-option
                 #include "pEqn.H"
-                rhoeps=rho*voidfraction;
             }
 
             if (pimple.turbCorr())
diff --git a/applications/solvers/cfdemSolverRhoPimple/pEqn.H b/applications/solvers/cfdemSolverRhoPimple/pEqn.H
index e107c5a4..d1c99774 100644
--- a/applications/solvers/cfdemSolverRhoPimple/pEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimple/pEqn.H
@@ -2,6 +2,7 @@ rho = thermo.rho();
 rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
+rhoeps = rho*voidfraction;
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
@@ -86,6 +87,8 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+rhoeps = rho*voidfraction;
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(voidfraction,p);