From 02a3327998728377ad75bb53e97366c2a9c460f3 Mon Sep 17 00:00:00 2001 From: andy Date: Wed, 8 May 2013 17:19:09 +0100 Subject: [PATCH] ENH: Extended forces and forceCoeffs to include porosity contributions --- .../forces/forceCoeffs/forceCoeffs.C | 4 +- .../functionObjects/forces/forces/forces.C | 219 ++++++++++++++---- .../functionObjects/forces/forces/forces.H | 17 +- 3 files changed, 190 insertions(+), 50 deletions(-) diff --git a/src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C b/src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C index 5a6d8ed4e7..9537f6a28a 100644 --- a/src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C +++ b/src/postProcessing/functionObjects/forces/forceCoeffs/forceCoeffs.C @@ -120,8 +120,8 @@ void Foam::forceCoeffs::write() scalar pDyn = 0.5*rhoRef_*magUInf_*magUInf_; - Field totForce(force_[0] + force_[1]); - Field totMoment(moment_[0] + moment_[1]); + Field totForce(force_[0] + force_[1] + force_[2]); + Field totMoment(moment_[0] + moment_[1] + moment_[2]); List > coeffs(3); coeffs[0].setSize(nBin_); diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C index fde287b8a4..c159b652e8 100644 --- a/src/postProcessing/functionObjects/forces/forces/forces.C +++ b/src/postProcessing/functionObjects/forces/forces/forces.C @@ -37,6 +37,7 @@ License #include "compressible/RAS/RASModel/RASModel.H" #include "compressible/LES/LESModel/LESModel.H" +#include "porosityModel.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -52,14 +53,15 @@ void Foam::forces::writeFileHeader(const label i) { file() << "# Time" << tab - << "forces(pressure, viscous) moment(pressure, viscous)"; + << "forces(pressure, viscous, porous) " + << "moment(pressure, viscous, porous)"; if (localSystem_) { file() << tab - << "local forces(pressure, viscous) " - << "local moment(pressure, viscous)"; + << "local forces(pressure, viscous, porous) " + << "local moment(pressure, viscous, porous)"; } file()<< endl; @@ -132,7 +134,7 @@ Foam::tmp Foam::forces::devRhoReff() const else { FatalErrorIn("forces::devRhoReff()") - << "No valid model for viscous stress calculation." + << "No valid model for viscous stress calculation" << exit(FatalError); return volSymmTensorField::null(); @@ -140,6 +142,46 @@ Foam::tmp Foam::forces::devRhoReff() const } +Foam::tmp Foam::forces::mu() const +{ + if (obr_.foundObject("thermophysicalProperties")) + { + const fluidThermo& thermo = + obr_.lookupObject("thermophysicalProperties"); + + return thermo.mu(); + } + else if + ( + obr_.foundObject("transportProperties") + ) + { + const singlePhaseTransportModel& laminarT = + obr_.lookupObject + ("transportProperties"); + + return rho()*laminarT.nu(); + } + else if (obr_.foundObject("transportProperties")) + { + const dictionary& transportProperties = + obr_.lookupObject("transportProperties"); + + dimensionedScalar nu(transportProperties.lookup("nu")); + + return rho()*nu; + } + else + { + FatalErrorIn("forces::mu()") + << "No valid model for dynamic viscosity calculation" + << exit(FatalError); + + return volScalarField::null(); + } +} + + Foam::tmp Foam::forces::rho() const { if (rhoName_ == "rhoInf") @@ -190,32 +232,35 @@ Foam::scalar Foam::forces::rho(const volScalarField& p) const void Foam::forces::applyBins ( - const label patchI, - const vectorField& fN, const vectorField& Md, - const vectorField& fT + const vectorField& fN, + const vectorField& fT, + const vectorField& fP, + const vectorField& d ) { if (nBin_ == 1) { force_[0][0] += sum(fN); force_[1][0] += sum(fT); - moment_[0][0] += sum(Md ^ fN); - moment_[1][0] += sum(Md ^ fT); + force_[2][0] += sum(fP); + moment_[0][0] += sum(Md^fN); + moment_[1][0] += sum(Md^fT); + moment_[2][0] += sum(Md^fP); } else { - const fvMesh& mesh = refCast(obr_); - const vectorField& Cf = mesh.C().boundaryField()[patchI]; - scalarField d((Cf & binDir_) - binMin_); + scalarField dd((d & binDir_) - binMin_); - forAll(d, i) + forAll(dd, i) { - label binI = floor(d[i]/binDx_); + label binI = floor(dd[i]/binDx_); force_[0][binI] += fN[i]; force_[1][binI] += fT[i]; + force_[2][binI] += fP[i]; moment_[0][binI] += Md[i]^fN[i]; moment_[1][binI] += Md[i]^fT[i]; + moment_[2][binI] += Md[i]^fP[i]; } } } @@ -239,7 +284,7 @@ void Foam::forces::writeBins() const Info<< " Writing bins to " << forcesDir << endl; } - wordList fieldNames(IStringStream("(pressure viscous)")()); + wordList fieldNames(IStringStream("(pressure viscous porous)")()); List > f(force_); List > m(moment_); @@ -250,8 +295,11 @@ void Foam::forces::writeBins() const { f[0][i] += f[0][i-1]; f[1][i] += f[1][i-1]; + f[2][i] += f[2][i-1]; + m[0][i] += m[0][i-1]; m[1][i] += m[1][i-1]; + m[2][i] += m[2][i-1]; } } @@ -264,12 +312,14 @@ void Foam::forces::writeBins() const if (localSystem_) { - List > localForce(2); - List > localMoment(2); + List > localForce(3); + List > localMoment(3); localForce[0] = coordSys_.localVector(force_[0]); localForce[1] = coordSys_.localVector(force_[1]); + localForce[2] = coordSys_.localVector(force_[2]); localMoment[0] = coordSys_.localVector(moment_[0]); localMoment[1] = coordSys_.localVector(moment_[1]); + localMoment[2] = coordSys_.localVector(moment_[2]); if (binCumulative_) { @@ -277,8 +327,10 @@ void Foam::forces::writeBins() const { localForce[0][i] += localForce[0][i-1]; localForce[1][i] += localForce[1][i-1]; + localForce[2][i] += localForce[2][i-1]; localMoment[0][i] += localMoment[0][i-1]; localMoment[1][i] += localMoment[1][i-1]; + localMoment[2][i] += localMoment[2][i-1]; } } @@ -306,8 +358,8 @@ Foam::forces::forces obr_(obr), active_(true), log_(false), - force_(2), - moment_(2), + force_(3), + moment_(3), patchSet_(), pName_(word::null), UName_(word::null), @@ -318,6 +370,7 @@ Foam::forces::forces pRef_(0), coordSys_(), localSystem_(false), + porosity_(false), nBin_(1), binDir_(vector::zero), binDx_(0.0), @@ -365,8 +418,8 @@ Foam::forces::forces obr_(obr), active_(true), log_(false), - force_(2), - moment_(2), + force_(3), + moment_(3), patchSet_(patchSet), pName_(pName), UName_(UName), @@ -377,6 +430,7 @@ Foam::forces::forces pRef_(pRef), coordSys_(coordSys), localSystem_(false), + porosity_(false), nBin_(1), binDir_(vector::zero), binDx_(0.0), @@ -475,6 +529,18 @@ void Foam::forces::read(const dictionary& dict) localSystem_ = true; } + if (dict.readIfPresent("porosity", porosity_)) + { + if (porosity_) + { + Info<< "Including porosity effects" << endl; + } + else + { + Info<< "Not including porosity effects" << endl; + } + } + if (dict.found("binData")) { const dictionary& binDict(dict.subDict("binData")); @@ -491,10 +557,11 @@ void Foam::forces::read(const dictionary& dict) else if ((nBin_ == 0) || (nBin_ == 1)) { nBin_ = 1; - force_[0].setSize(1); - force_[1].setSize(1); - moment_[0].setSize(1); - moment_[1].setSize(1); + forAll(force_, i) + { + force_[i].setSize(1); + moment_[i].setSize(1); + } } if (nBin_ > 1) @@ -546,8 +613,10 @@ void Foam::forces::read(const dictionary& dict) // allocate storage for forces and moments force_[0].setSize(1); force_[1].setSize(1); + force_[2].setSize(1); moment_[0].setSize(1); moment_[1].setSize(1); + moment_[2].setSize(1); } } } @@ -581,28 +650,50 @@ void Foam::forces::write() if (log_) { Info<< type() << " output:" << nl - << " forces(pressure,viscous)" - << "(" << sum(force_[0]) << "," << sum(force_[1]) << ")" << nl - << " moment(pressure,viscous)" - << "(" << sum(moment_[0]) << "," << sum(moment_[1]) << ")" + << " forces(pressure,viscous,porous) = (" + << sum(force_[0]) << "," + << sum(force_[1]) << "," + << sum(force_[2]) << ")" << nl + << " moment(pressure,viscous,porous) = (" + << sum(moment_[0]) << "," + << sum(moment_[1]) << "," + << sum(moment_[2]) << ")" << nl; } file() << obr_.time().value() << tab - << "(" << sum(force_[0]) << "," << sum(force_[1]) << ") " - << "(" << sum(moment_[0]) << "," << sum(moment_[1]) << ")" + << "(" + << sum(force_[0]) << "," + << sum(force_[1]) << "," + << sum(force_[2]) + << ") " + << "(" + << sum(moment_[0]) << "," + << sum(moment_[1]) << "," + << sum(moment_[2]) + << ")" << endl; if (localSystem_) { - vectorField localForceP(coordSys_.localVector(force_[0])); - vectorField localForceV(coordSys_.localVector(force_[1])); - vectorField localMomentP(coordSys_.localVector(moment_[0])); - vectorField localMomentV(coordSys_.localVector(moment_[1])); + vectorField localForceN(coordSys_.localVector(force_[0])); + vectorField localForceT(coordSys_.localVector(force_[1])); + vectorField localForceP(coordSys_.localVector(force_[2])); + vectorField localMomentN(coordSys_.localVector(moment_[0])); + vectorField localMomentT(coordSys_.localVector(moment_[1])); + vectorField localMomentP(coordSys_.localVector(moment_[2])); file() << obr_.time().value() << tab - << "(" << sum(localForceP) << "," << sum(localForceV) << ") " - << "(" << sum(localMomentP) << "," << sum(localMomentV) << ")" + << "(" + << sum(localForceN) << "," + << sum(localForceT) << "," + << sum(localForceP) + << ") " + << "(" + << sum(localMomentN) << "," + << sum(localMomentT) << "," + << sum(localMomentP) + << ")" << endl; } @@ -620,9 +711,11 @@ void Foam::forces::calcForcesMoment() { force_[0] = vector::zero; force_[1] = vector::zero; + force_[2] = vector::zero; moment_[0] = vector::zero; moment_[1] = vector::zero; + moment_[2] = vector::zero; if (directForceDensity_) { @@ -656,7 +749,10 @@ void Foam::forces::calcForcesMoment() // Tangential force (total force minus normal fN) vectorField fT(sA*fD.boundaryField()[patchI] - fN); - applyBins(patchI, fN, Md, fT); + //- Porous force + vectorField fP(Md.size(), vector::zero); + + applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]); } } else @@ -685,14 +781,51 @@ void Foam::forces::calcForcesMoment() mesh.C().boundaryField()[patchI] - coordSys_.origin() ); - vectorField pf + vectorField fN ( rho(p)*Sfb[patchI]*(p.boundaryField()[patchI] - pRef) ); - vectorField vf(Sfb[patchI] & devRhoReffb[patchI]); + vectorField fT(Sfb[patchI] & devRhoReffb[patchI]); - applyBins(patchI, pf, Md, vf); + vectorField fP(Md.size(), vector::zero); + + applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]); + } + } + + if (porosity_) + { + const volVectorField& U = obr_.lookupObject(UName_); + const volScalarField rho(this->rho()); + const volScalarField mu(this->mu()); + + const fvMesh& mesh = U.mesh(); + + const HashTable models = + obr_.lookupClass(); + + forAllConstIter(HashTable, models, iter) + { + const porosityModel& pm = *iter(); + + vectorField fPTot(pm.force(U, rho, mu)); + + const labelList& cellZoneIDs = pm.cellZoneIDs(); + + forAll(cellZoneIDs, i) + { + label zoneI = cellZoneIDs[i]; + const cellZone& cZone = mesh.cellZones()[zoneI]; + + const vectorField d(mesh.C(), cZone); + const vectorField fP(fPTot, cZone); + const vectorField Md(d - coordSys_.origin()); + + const vectorField fDummy(Md.size(), vector::zero); + + applyBins(Md, fDummy, fDummy, fP, d); + } } } @@ -703,13 +836,13 @@ void Foam::forces::calcForcesMoment() Foam::vector Foam::forces::forceEff() const { - return sum(force_[0]) + sum(force_[1]); + return sum(force_[0]) + sum(force_[1]) + sum(force_[2]); } Foam::vector Foam::forces::momentEff() const { - return sum(moment_[0]) + sum(moment_[1]); + return sum(moment_[0]) + sum(moment_[1]) + sum(moment_[2]); } diff --git a/src/postProcessing/functionObjects/forces/forces/forces.H b/src/postProcessing/functionObjects/forces/forces/forces.H index c0f60e614a..dd87184c2f 100644 --- a/src/postProcessing/functionObjects/forces/forces/forces.H +++ b/src/postProcessing/functionObjects/forces/forces/forces.H @@ -149,10 +149,10 @@ protected: //- Switch to send output to Info as well as to file Switch log_; - //- Pressure and viscous force per bin + //- Pressure, viscous and porous force per bin List > force_; - //- Pressure and viscous pressure per bin + //- Pressure, viscous and porous moment per bin List > moment_; @@ -188,6 +188,9 @@ protected: //- Flag to indicate whether we are using a local co-ordinate sys bool localSystem_; + //- Flag to include porosity effects + bool porosity_; + // Bin information @@ -221,6 +224,9 @@ protected: //- Return the effective viscous stress (laminar + turbulent). tmp devRhoReff() const; + //- Dynamic viscosity field + tmp mu() const; + //- Return rho if rhoName is specified otherwise rhoRef tmp rho() const; @@ -231,10 +237,11 @@ protected: //- Accumulate bin data void applyBins ( - const label patchI, - const vectorField& fN, const vectorField& Md, - const vectorField& fT + const vectorField& fN, + const vectorField& fT, + const vectorField& fP, + const vectorField& d ); //- Helper function to write bin data