diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C index 462116ca4a..97f68bb216 100644 --- a/src/postProcessing/functionObjects/forces/forces/forces.C +++ b/src/postProcessing/functionObjects/forces/forces/forces.C @@ -190,6 +190,7 @@ Foam::forces::forces directForceDensity_(false), fDName_(""), rhoRef_(VGREAT), + pRef_(0), CofR_(vector::zero), forcesFilePtr_(NULL) { @@ -282,13 +283,15 @@ void Foam::forces::read(const dictionary& dict) Info<< " or " << rhoName_; } - Info<< " in database." << nl - << " De-activating forces." + Info<< " in database." << nl << " De-activating forces." << endl; } // Reference density needed for incompressible calculations rhoRef_ = readScalar(dict.lookup("rhoInf")); + + // Reference pressure, 0 by default + pRef_ = dict.lookupOrDefault("pRef", 0.0); } // Centre of rotation for moment calculations @@ -443,13 +446,16 @@ Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const const volSymmTensorField::GeometricBoundaryField& devRhoReffb = tdevRhoReff().boundaryField(); + // Scale pRef by density for incompressible simulations + scalar pRef = pRef_/rho(p); + forAllConstIter(labelHashSet, patchSet_, iter) { label patchi = iter.key(); vectorField Md = mesh.C().boundaryField()[patchi] - CofR_; - vectorField pf = Sfb[patchi]*p.boundaryField()[patchi]; + vectorField pf = Sfb[patchi]*(p.boundaryField()[patchi] - pRef); fm.first().first() += rho(p)*sum(pf); fm.second().first() += rho(p)*sum(Md ^ pf); diff --git a/src/postProcessing/functionObjects/forces/forces/forces.H b/src/postProcessing/functionObjects/forces/forces/forces.H index 6fcdccc50e..42bb67c23b 100644 --- a/src/postProcessing/functionObjects/forces/forces/forces.H +++ b/src/postProcessing/functionObjects/forces/forces/forces.H @@ -144,6 +144,9 @@ protected: //- Reference density needed for incompressible calculations scalar rhoRef_; + //- Reference pressure + scalar pRef_; + //- Centre of rotation vector CofR_;