diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C
index e779aac4..73ec3fca 100644
--- a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C
+++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C
@@ -38,6 +38,7 @@ Description
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
+#include "fvOptions.H"
#include "cfdemCloud.H"
#include "implicitCouple.H"
@@ -54,6 +55,7 @@ int main(int argc, char *argv[])
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
+ #include "createFvOptions.H"
#include "initContinuityErrs.H"
// create cfdemCloud
@@ -101,78 +103,19 @@ int main(int argc, char *argv[])
// Pressure-velocity PISO corrector
{
// Momentum predictor
- fvVectorMatrix UEqn
- (
- fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
- + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
-// + turbulence->divDevReff(U)
- + particleCloud.divVoidfractionTau(U, voidfraction)
- ==
- - fvm::Sp(Ksl/rho,U)
- );
-
- UEqn.relax();
- if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
- solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
- else if (piso.momentumPredictor())
- solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+ #include "UEqn.H"
// --- PISO loop
while (piso.correct())
{
- volScalarField rUA = 1.0/UEqn.A();
-
- surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
- volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
- surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
-
- U = rUA*UEqn.H();
-
-
- phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
- + rUAfvoidfraction*fvc::ddtCorr(U, phi);
-
- surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
- surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
- if (modelType=="A")
- rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
- // Non-orthogonal pressure corrector loop
- while (piso.correctNonOrthogonal())
- {
- // Pressure corrector
- fvScalarMatrix pEqn
- (
- fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
- );
- pEqn.setReference(pRefCell, pRefValue);
-
- pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
-
- if (piso.finalNonOrthogonalIter())
- {
- phiGes -= pEqn.flux();
- phi = phiGes;
- }
-
- } // end non-orthogonal corrector loop
-
- #include "continuityErrorPhiPU.H"
-
- if (modelType=="B" || modelType=="Bfull")
- U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
- else
- U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
- U.correctBoundaryConditions();
-
- } // end piso loop
+ #include "pEqn.H"
+ }
}
+ laminarTransport.correct();
turbulence->correct();
- }// end solveFlow
+ }
else
{
Info << "skipping flow solution." << endl;
diff --git a/applications/solvers/cfdemSolverPiso/createFields.H b/applications/solvers/cfdemSolverPiso/createFields.H
index 849ef338..8f0b5584 100644
--- a/applications/solvers/cfdemSolverPiso/createFields.H
+++ b/applications/solvers/cfdemSolverPiso/createFields.H
@@ -122,3 +122,5 @@ surfaceScalarField phi
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
+
+#include "createMRF.H"
\ No newline at end of file
diff --git a/applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C b/applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C
index 0d8e5137..c8322fad 100644
--- a/applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C
+++ b/applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C
@@ -25,12 +25,12 @@ License
along with CFDEMcoupling. If not, see .
Application
- cfdemSolverPisoMS
+ cfdemSolverPiso
Description
Transient solver for incompressible flow.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
- The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
+ The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
where additional functionality for CFD-DEM coupling is added.
\*---------------------------------------------------------------------------*/
@@ -38,11 +38,13 @@ Description
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
+#include "fvOptions.H"
#include "cfdemCloudMS.H"
#include "implicitCouple.H"
#include "clockModel.H"
#include "smoothingModel.H"
+#include "forceModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -53,6 +55,7 @@ int main(int argc, char *argv[])
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
+ #include "createFvOptions.H"
#include "initContinuityErrs.H"
// create cfdemCloud
@@ -61,11 +64,10 @@ int main(int argc, char *argv[])
#include "checkModelType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+ Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
- Info<< "\nStarting time loop\n" << endl;
- particleCloud.clockM().start(1,"Global");
+ particleCloud.clockM().start(1,"Global");
Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -73,91 +75,51 @@ int main(int argc, char *argv[])
// do particle stuff
particleCloud.clockM().start(2,"Coupling");
- particleCloud.evolve(voidfraction,Us,U);
-
+ bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+ if(hasEvolved)
+ {
+ particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+ }
+
Info << "update Ksl.internalField()" << endl;
Ksl = particleCloud.momCoupleM(0).impMomSource();
- particleCloud.smoothingM().smoothen(Ksl);
Ksl.correctBoundaryConditions();
+ //Force Checks
+ vector fTotal(0,0,0);
+ vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+ reduce(fImpTotal, sumOp());
+ Info << "TotalForceExp: " << fTotal << endl;
+ Info << "TotalForceImp: " << fImpTotal << endl;
+
#include "solverDebugInfo.H"
particleCloud.clockM().stop("Coupling");
particleCloud.clockM().start(26,"Flow");
- // Pressure-velocity PISO corrector
+
+ if(particleCloud.solveFlow())
{
- // Momentum predictor
- fvVectorMatrix UEqn
- (
- fvm::ddt(voidfraction,U) //particleCloud.ddtVoidfractionU(U,voidfraction) //
- + fvm::div(phi, U)
-// + turbulence->divDevReff(U)
- + particleCloud.divVoidfractionTau(U, voidfraction)
- ==
- - fvm::Sp(Ksl/rho,U)
- );
-
- if (modelType=="B")
- UEqn == - fvc::grad(p) + Ksl/rho*Us;
- else
- UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us;
-
- UEqn.relax();
-
- if (piso.momentumPredictor())
- solve(UEqn);
-
- // --- PISO loop
-
- while (piso.correct())
+ // Pressure-velocity PISO corrector
{
- volScalarField rUA = 1.0/UEqn.A();
+ // Momentum predictor
+ #include "UEqn.H"
- surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
- volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+ // --- PISO loop
- U = rUA*UEqn.H();
-
- phi = (fvc::interpolate(U*voidfraction) & mesh.Sf() );
- //+ fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
- surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
- surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
- if (modelType=="A")
- rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
- // Non-orthogonal pressure corrector loop
- while (piso.correctNonOrthogonal())
+ while (piso.correct())
{
- // Pressure corrector
- fvScalarMatrix pEqn
- (
- fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
- );
- pEqn.setReference(pRefCell, pRefValue);
+ #include "pEqn.H"
+ }
+ }
- pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
-
- if (piso.finalNonOrthogonalIter())
- {
- phiGes -= pEqn.flux();
- }
-
- } // end non-orthogonal corrector loop
-
- #include "continuityErrorPhiPU.H"
-
- if (modelType=="B")
- U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
- else
- U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
- U.correctBoundaryConditions();
-
- } // end piso loop
+ laminarTransport.correct();
+ turbulence->correct();
+ }
+ else
+ {
+ Info << "skipping flow solution." << endl;
}
-
- turbulence->correct();
runTime.write();
@@ -170,7 +132,7 @@ int main(int argc, char *argv[])
}
Info<< "End\n" << endl;
-
+
return 0;
}
diff --git a/applications/solvers/cfdemSolverPisoScalar/Make/options b/applications/solvers/cfdemSolverPisoScalar/Make/options
index 91f9f553..69598800 100644
--- a/applications/solvers/cfdemSolverPisoScalar/Make/options
+++ b/applications/solvers/cfdemSolverPisoScalar/Make/options
@@ -7,6 +7,7 @@ EXE_INC = \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+ -I ../cfdemSolverPiso \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
diff --git a/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C b/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C
index d2621d91..788c8f02 100644
--- a/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C
+++ b/applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C
@@ -25,7 +25,7 @@ License
along with CFDEMcoupling. If not, see .
Application
- cfdemSolverPisoScalar
+ cfdemSolverPiso
Description
Transient solver for incompressible flow.
@@ -38,9 +38,11 @@ Description
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
+#include "fvOptions.H"
#include "cfdemCloud.H"
#include "implicitCouple.H"
+#include "clockModel.H"
#include "smoothingModel.H"
#include "forceModel.H"
@@ -53,6 +55,7 @@ int main(int argc, char *argv[])
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
+ #include "createFvOptions.H"
#include "initContinuityErrs.H"
// create cfdemCloud
@@ -64,11 +67,14 @@ int main(int argc, char *argv[])
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
+ particleCloud.clockM().start(1,"Global");
+
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "CourantNo.H"
// do particle stuff
+ particleCloud.clockM().start(2,"Coupling");
bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
if(hasEvolved)
@@ -80,100 +86,38 @@ int main(int argc, char *argv[])
Ksl = particleCloud.momCoupleM(0).impMomSource();
Ksl.correctBoundaryConditions();
+ //Force Checks
+ vector fTotal(0,0,0);
+ vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+ reduce(fImpTotal, sumOp());
+ Info << "TotalForceExp: " << fTotal << endl;
+ Info << "TotalForceImp: " << fImpTotal << endl;
#include "solverDebugInfo.H"
+ particleCloud.clockM().stop("Coupling");
- // get scalar source from DEM
- particleCloud.forceM(1).manipulateScalarField(Tsource);
- Tsource.correctBoundaryConditions();
-
- // solve scalar transport equation
- fvScalarMatrix TEqn
- (
- fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
- + fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
- - fvm::laplacian(DT*voidfraction, T)
- ==
- Tsource
- );
- TEqn.relax();
- TEqn.solve();
+ particleCloud.clockM().start(26,"Flow");
+
+ #include "TEqn.H"
if(particleCloud.solveFlow())
{
// Pressure-velocity PISO corrector
{
// Momentum predictor
- fvVectorMatrix UEqn
- (
- fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
- + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
-// + turbulence->divDevReff(U)
- + particleCloud.divVoidfractionTau(U, voidfraction)
- ==
- - fvm::Sp(Ksl/rho,U)
- );
-
- UEqn.relax();
- if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
- solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
- else if (piso.momentumPredictor())
- solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+ #include "UEqn.H"
// --- PISO loop
while (piso.correct())
{
- volScalarField rUA = 1.0/UEqn.A();
-
- surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
- volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
- surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
-
- U = rUA*UEqn.H();
-
- phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
- + rUAfvoidfraction*fvc::ddtCorr(U, phi);
-
- surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
- surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
- if (modelType=="A")
- rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
- while (piso.correctNonOrthogonal())
- {
- // Pressure corrector
- fvScalarMatrix pEqn
- (
- fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
- );
- pEqn.setReference(pRefCell, pRefValue);
-
- pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
-
- if (piso.finalNonOrthogonalIter())
- {
- phiGes -= pEqn.flux();
- phi = phiGes;
- }
-
- } // end non-orthogonal corrector loop
-
- #include "continuityErrorPhiPU.H"
-
- if (modelType=="B" || modelType=="Bfull")
- U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
- else
- U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
- U.correctBoundaryConditions();
-
- } // end piso loop
+ #include "pEqn.H"
+ }
}
+ laminarTransport.correct();
turbulence->correct();
- }// end solveFlow
+ }
else
{
Info << "skipping flow solution." << endl;
@@ -184,6 +128,9 @@ int main(int argc, char *argv[])
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
+
+ particleCloud.clockM().stop("Flow");
+ particleCloud.clockM().stop("Global");
}
Info<< "End\n" << endl;
diff --git a/applications/solvers/cfdemSolverPisoScalar/createFields.H b/applications/solvers/cfdemSolverPisoScalar/createFields.H
index bff7514c..fe5dd29a 100644
--- a/applications/solvers/cfdemSolverPisoScalar/createFields.H
+++ b/applications/solvers/cfdemSolverPisoScalar/createFields.H
@@ -1,36 +1,36 @@
- Info<< "Reading field p\n" << endl;
- volScalarField p
- (
- IOobject
- (
- "p",
- runTime.timeName(),
- mesh,
- IOobject::MUST_READ,
- IOobject::AUTO_WRITE
- ),
- mesh
- );
-
- Info<< "Reading physical velocity field U" << endl;
- Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
- volVectorField U
- (
- IOobject
- (
- "U",
- runTime.timeName(),
- mesh,
- IOobject::MUST_READ,
- IOobject::AUTO_WRITE
- ),
- mesh
- );
-
-//========================
-// drag law modelling
-//========================
-
+ Info<< "Reading field p\n" << endl;
+ volScalarField p
+ (
+ IOobject
+ (
+ "p",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ );
+
+ Info<< "Reading physical velocity field U" << endl;
+ Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
+ volVectorField U
+ (
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ );
+
+//========================
+// drag law modelling
+//========================
+
Info<< "\nReading momentum exchange field Ksl\n" << endl;
volScalarField Ksl
(
@@ -44,8 +44,8 @@
),
mesh
//dimensionedScalar("0", dimensionSet(0, 0, -1, 0, 0), 1.0)
- );
-
+ );
+
Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
volScalarField voidfraction
(
@@ -58,8 +58,8 @@
IOobject::AUTO_WRITE
),
mesh
- );
-
+ );
+
Info<< "\nCreating density field rho\n" << endl;
volScalarField rho
(
@@ -71,27 +71,27 @@
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
- mesh,
+ mesh,
dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
- );
-
- Info<< "Reading particle velocity field Us\n" << endl;
- volVectorField Us
- (
- IOobject
- (
- "Us",
- runTime.timeName(),
- mesh,
- IOobject::MUST_READ,
- IOobject::AUTO_WRITE
- ),
- mesh
- );
-
-//========================
-// scalar field modelling
-//========================
+ );
+
+ Info<< "Reading particle velocity field Us\n" << endl;
+ volVectorField Us
+ (
+ IOobject
+ (
+ "Us",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ );
+
+//========================
+// scalar field modelling
+//========================
Info<< "\nCreating dummy density field rho = 1\n" << endl;
volScalarField T
(
@@ -103,10 +103,10 @@
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
- mesh//,
+ mesh//,
//dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 273.15)
- );
-
+ );
+
Info<< "\nCreating fluid-particle heat flux field\n" << endl;
volScalarField Tsource
(
@@ -118,57 +118,59 @@
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
- mesh//,
+ mesh//,
//dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 0.0)
- );
-
- IOdictionary transportProperties
- (
- IOobject
- (
- "transportProperties",
- runTime.constant(),
- mesh,
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- )
- );
-
- dimensionedScalar DT
- (
- transportProperties.lookup("DT")
- );
-
-//========================
-
-//# include "createPhi.H"
-#ifndef createPhi_H
-#define createPhi_H
-Info<< "Reading/calculating face flux field phi\n" << endl;
-surfaceScalarField phi
- (
- IOobject
- (
- "phi",
- runTime.timeName(),
- mesh,
- IOobject::READ_IF_PRESENT,
- IOobject::AUTO_WRITE
- ),
- linearInterpolate(U*voidfraction) & mesh.Sf()
- );
-#endif
-
-
-
- label pRefCell = 0;
- scalar pRefValue = 0.0;
- setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
- singlePhaseTransportModel laminarTransport(U, phi);
-
- autoPtr turbulence
- (
- incompressible::turbulenceModel::New(U, phi, laminarTransport)
- );
+ );
+
+ IOdictionary transportProperties
+ (
+ IOobject
+ (
+ "transportProperties",
+ runTime.constant(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
+ );
+
+ dimensionedScalar DT
+ (
+ transportProperties.lookup("DT")
+ );
+
+//========================
+
+//# include "createPhi.H"
+#ifndef createPhi_H
+#define createPhi_H
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+ (
+ IOobject
+ (
+ "phi",
+ runTime.timeName(),
+ mesh,
+ IOobject::READ_IF_PRESENT,
+ IOobject::AUTO_WRITE
+ ),
+ linearInterpolate(U*voidfraction) & mesh.Sf()
+ );
+#endif
+
+
+
+ label pRefCell = 0;
+ scalar pRefValue = 0.0;
+ setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+ singlePhaseTransportModel laminarTransport(U, phi);
+
+ autoPtr turbulence
+ (
+ incompressible::turbulenceModel::New(U, phi, laminarTransport)
+ );
+
+#include "createMRF.H"
\ No newline at end of file
diff --git a/src/lagrangian/cfdemParticle/cfdTools/continuityErrorPhiPU.H b/src/lagrangian/cfdemParticle/cfdTools/continuityErrorPhiPU.H
index c58afcfb..6506ef64 100644
--- a/src/lagrangian/cfdemParticle/cfdTools/continuityErrorPhiPU.H
+++ b/src/lagrangian/cfdemParticle/cfdTools/continuityErrorPhiPU.H
@@ -34,7 +34,7 @@ Description
\*---------------------------------------------------------------------------*/
{
- volScalarField contErr( fvc::div(phiGes) + fvc::ddt(voidfraction) );
+ volScalarField contErr( fvc::div(phi) + fvc::ddt(voidfraction) );
scalar sumLocalContErr = runTime.deltaTValue()*
mag(contErr)().weightedAverage(mesh.V()).value();