From 3d2d26334a62c99ce17498541357a8f22a391320 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Wed, 11 Jan 2017 11:57:48 +0000 Subject: [PATCH] Bugfix: added ramping to MRF zones --- .../rhoPorousMRFPimpleFoam/pEqn.H | 3 + .../steadyCompressibleMRFFoam/pEqn.H | 3 + .../steadyUniversalMRFFoam/pEqn.H | 3 + .../incompressible/MRFSimpleFoam/pEqn.H | 4 + .../solvers/multiphase/MRFInterFoam/pEqn.H | 3 + .../cfdTools/general/MRF/MRFZone.C | 116 ++++++++++-------- .../cfdTools/general/MRF/MRFZones.C | 26 ++-- .../cfdTools/general/porousMedia/porousZone.H | 17 ++- .../general/porousMedia/porousZones.H | 2 +- 9 files changed, 104 insertions(+), 73 deletions(-) diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H index 920ba0e5d..18ee1a2b4 100644 --- a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H @@ -3,6 +3,9 @@ rho = thermo.rho(); volScalarField rUA = 1.0/UEqn().A(); U = rUA*UEqn().H(); +// Update boundary velocity for consistency with the flux +mrfZones.correctBoundaryVelocity(U); + if (pimple.nCorrPISO() <= 1) { UEqn.clear(); diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H index ff35bebaa..61fed64e8 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H @@ -12,6 +12,9 @@ { U = rUA*UEqn.H(); + // Update boundary velocity for consistency with the flux + mrfZones.correctBoundaryVelocity(U); + // Calculate phi for boundary conditions phi = rhof*fvc::interpolate(U) & mesh.Sf(); diff --git a/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H b/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H index 1d4b2dc94..bba75ac7f 100644 --- a/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyUniversalMRFFoam/pEqn.H @@ -12,6 +12,9 @@ { U = rUA*UEqn.H(); + // Update boundary velocity for consistency with the flux + mrfZones.correctBoundaryVelocity(U); + // Calculate phi for boundary conditions phi = rhof*fvc::interpolate(U) & mesh.Sf(); diff --git a/applications/solvers/incompressible/MRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/MRFSimpleFoam/pEqn.H index 1752fc663..d5bb1e361 100644 --- a/applications/solvers/incompressible/MRFSimpleFoam/pEqn.H +++ b/applications/solvers/incompressible/MRFSimpleFoam/pEqn.H @@ -9,6 +9,10 @@ 1/HUEqn().A() ); + // Update boundary velocity for consistency with the flux + // This is needed for a ramped MRF + mrfZones.correctBoundaryVelocity(U); + // Store velocity under-relaxation point before using U for // the flux precursor U.storePrevIter(); diff --git a/applications/solvers/multiphase/MRFInterFoam/pEqn.H b/applications/solvers/multiphase/MRFInterFoam/pEqn.H index 49d6dc579..eebfdfa7d 100644 --- a/applications/solvers/multiphase/MRFInterFoam/pEqn.H +++ b/applications/solvers/multiphase/MRFInterFoam/pEqn.H @@ -4,6 +4,9 @@ U = rUA*UEqn.H(); + // Update boundary velocity for consistency with the flux + mrfZones.correctBoundaryVelocity(U); + surfaceScalarField phiU ( "phiU", diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C index 09dfb6e25..a26ff72ce 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C @@ -32,6 +32,8 @@ License #include "faceSet.H" #include "geometricOneField.H" +using namespace Foam::mathematicalConstant; + // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // defineTypeNameAndDebug(Foam::MRFZone, 0); @@ -44,9 +46,9 @@ void Foam::MRFZone::setMRFFaces() const polyBoundaryMesh& patches = mesh_.boundaryMesh(); // Type per face: - // 0:not in zone - // 1:moving with frame - // 2:other + // 0: not in zone + // 1: moving with frame + // 2: other labelList faceType(mesh_.nFaces(), 0); // Determine faces in cell zone @@ -68,7 +70,6 @@ void Foam::MRFZone::setMRFFaces() } } - label nZoneFaces = 0; for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) @@ -83,7 +84,7 @@ void Foam::MRFZone::setMRFFaces() labelHashSet excludedPatches(excludedPatchLabels_); - + Info<< "Excluded patches: " << excludedPatchLabels_ << endl; forAll (patches, patchI) { const polyPatch& pp = patches[patchI]; @@ -142,9 +143,9 @@ void Foam::MRFZone::setMRFFaces() { const polyPatch& pp = patches[patchi]; - forAll (pp, patchFacei) + forAll (pp, patchFaceI) { - label faceI = pp.start()+patchFacei; + label faceI = pp.start() + patchFaceI; if (faceType[faceI] == 1) { @@ -159,29 +160,29 @@ void Foam::MRFZone::setMRFFaces() includedFaces_.setSize(patches.size()); excludedFaces_.setSize(patches.size()); - forAll (nIncludedFaces, patchi) + forAll (nIncludedFaces, patchI) { - includedFaces_[patchi].setSize(nIncludedFaces[patchi]); - excludedFaces_[patchi].setSize(nExcludedFaces[patchi]); + includedFaces_[patchI].setSize(nIncludedFaces[patchI]); + excludedFaces_[patchI].setSize(nExcludedFaces[patchI]); } nIncludedFaces = 0; nExcludedFaces = 0; - forAll (patches, patchi) + forAll (patches, patchI) { - const polyPatch& pp = patches[patchi]; + const polyPatch& pp = patches[patchI]; - forAll (pp, patchFacei) + forAll (pp, patchFaceI) { - label faceI = pp.start()+patchFacei; + label faceI = pp.start() + patchFaceI; if (faceType[faceI] == 1) { - includedFaces_[patchi][nIncludedFaces[patchi]++] = patchFacei; + includedFaces_[patchI][nIncludedFaces[patchI]++] = patchFaceI; } else if (faceType[faceI] == 2) { - excludedFaces_[patchi][nExcludedFaces[patchi]++] = patchFacei; + excludedFaces_[patchI][nExcludedFaces[patchI]++] = patchFaceI; } } } @@ -201,12 +202,13 @@ void Foam::MRFZone::setMRFFaces() internalFaces.write(); faceSet MRFFaces(mesh_, "includedFaces", 100); - forAll (includedFaces_, patchi) + + forAll (includedFaces_, patchI) { - forAll (includedFaces_[patchi], i) + forAll (includedFaces_[patchI], i) { - label patchFacei = includedFaces_[patchi][i]; - MRFFaces.insert(patches[patchi].start()+patchFacei); + label patchFaceI = includedFaces_[patchI][i]; + MRFFaces.insert(patches[patchI].start() + patchFaceI); } } Pout<< "Writing " << MRFFaces.size() @@ -215,12 +217,12 @@ void Foam::MRFZone::setMRFFaces() MRFFaces.write(); faceSet excludedFaces(mesh_, "excludedFaces", 100); - forAll (excludedFaces_, patchi) + forAll (excludedFaces_, patchI) { - forAll (excludedFaces_[patchi], i) + forAll (excludedFaces_[patchI], i) { - label patchFacei = excludedFaces_[patchi][i]; - excludedFaces.insert(patches[patchi].start() + patchFacei); + label patchFaceI = excludedFaces_[patchI][i]; + excludedFaces.insert(patches[patchI].start() + patchFaceI); } } Pout<< "Writing " << excludedFaces.size() @@ -240,11 +242,15 @@ Foam::vector Foam::MRFZone::Omega() const else { // Ramping - Info<< "MRF Ramp: " << Foam::min(mesh_.time().value()/rampTime_, 1) - << endl; + const scalar t = mesh_.time().value(); + const scalar ramp = sin(2*pi/(4*rampTime_)*Foam::min(rampTime_, t)); - return Foam::min(mesh_.time().value()/rampTime_, 1)* - omega_.value()*axis_.value(); + Info<< "MRF Ramp: " << ramp << endl; + + return ramp*omega_.value()*axis_.value(); + + // return Foam::min(mesh_.time().value()/rampTime_, 1.0)* + // omega_.value()*axis_.value(); } } @@ -304,7 +310,6 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is) } } - bool cellZoneFound = (cellZoneID_ != -1); reduce(cellZoneFound, orOp()); @@ -317,6 +322,11 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is) << exit(FatalError); } + + Info<< "Creating MRF for cell zone " << name_ << ". rpm = " + << 30*omega_.value()/mathematicalConstant::pi + << endl; + setMRFFaces(); } @@ -384,24 +394,24 @@ void Foam::MRFZone::relativeVelocity(volVectorField& U) const U[celli] -= (rotVel ^ (C[celli] - origin)); } - // Included patches + // Included faces forAll (includedFaces_, patchi) { forAll (includedFaces_[patchi], i) { - label patchFacei = includedFaces_[patchi][i]; - U.boundaryField()[patchi][patchFacei] = vector::zero; + label patchFaceI = includedFaces_[patchi][i]; + U.boundaryField()[patchi][patchFaceI] = vector::zero; } } - // Excluded patches + // Excluded faces forAll (excludedFaces_, patchi) { forAll (excludedFaces_[patchi], i) { - label patchFacei = excludedFaces_[patchi][i]; - U.boundaryField()[patchi][patchFacei] -= - (rotVel ^ (C.boundaryField()[patchi][patchFacei] - origin)); + label patchFaceI = excludedFaces_[patchi][i]; + U.boundaryField()[patchi][patchFaceI] -= + (rotVel ^ (C.boundaryField()[patchi][patchFaceI] - origin)); } } } @@ -422,25 +432,25 @@ void Foam::MRFZone::absoluteVelocity(volVectorField& U) const U[celli] += (rotVel ^ (C[celli] - origin)); } - // Included patches + // Included faces forAll (includedFaces_, patchi) { forAll (includedFaces_[patchi], i) { - label patchFacei = includedFaces_[patchi][i]; - U.boundaryField()[patchi][patchFacei] = - (rotVel ^ (C.boundaryField()[patchi][patchFacei] - origin)); + label patchFaceI = includedFaces_[patchi][i]; + U.boundaryField()[patchi][patchFaceI] = + (rotVel ^ (C.boundaryField()[patchi][patchFaceI] - origin)); } } - // Excluded patches + // Excluded faces forAll (excludedFaces_, patchi) { forAll (excludedFaces_[patchi], i) { - label patchFacei = excludedFaces_[patchi][i]; - U.boundaryField()[patchi][patchFacei] += - (rotVel ^ (C.boundaryField()[patchi][patchFacei] - origin)); + label patchFaceI = excludedFaces_[patchi][i]; + U.boundaryField()[patchi][patchFaceI] += + (rotVel ^ (C.boundaryField()[patchi][patchFaceI] - origin)); } } } @@ -499,10 +509,10 @@ void Foam::MRFZone::meshPhi { forAll (includedFaces_[patchi], i) { - label patchFacei = includedFaces_[patchi][i]; + label patchFaceI = includedFaces_[patchi][i]; - phi.boundaryField()[patchi][patchFacei] = - (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)); + phi.boundaryField()[patchi][patchFaceI] = + (rotVel ^ (Cf.boundaryField()[patchi][patchFaceI] - origin)); } } @@ -511,10 +521,10 @@ void Foam::MRFZone::meshPhi { forAll (excludedFaces_[patchi], i) { - label patchFacei = excludedFaces_[patchi][i]; + label patchFaceI = excludedFaces_[patchi][i]; - phi.boundaryField()[patchi][patchFacei] = - (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)); + phi.boundaryField()[patchi][patchFaceI] = + (rotVel ^ (Cf.boundaryField()[patchi][patchFaceI] - origin)); } } } @@ -523,7 +533,7 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const { const vector& origin = origin_.value(); const vector rotVel = Omega(); - + Info<< "rotVel: " << rotVel << endl; // Included patches forAll (includedFaces_, patchi) { @@ -533,9 +543,9 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const forAll (includedFaces_[patchi], i) { - label patchFacei = includedFaces_[patchi][i]; + label patchFaceI = includedFaces_[patchi][i]; - pfld[patchFacei] = (rotVel ^ (patchC[patchFacei] - origin)); + pfld[patchFaceI] = (rotVel ^ (patchC[patchFaceI] - origin)); } U.boundaryField()[patchi] == pfld; diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C index b2a6a780c..e5f29c690 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C @@ -77,7 +77,7 @@ Foam::tmp Foam::MRFZones::fluxCorrection() const ); surfaceScalarField& MRFZonesPhiCorr = tMRFZonesPhiCorr(); - forAll(*this, i) + forAll (*this, i) { operator[](i).relativeFlux(MRFZonesPhiCorr); } @@ -106,7 +106,7 @@ Foam::tmp Foam::MRFZones::meshPhi() const ); surfaceVectorField& MRFZonesPhiCorr = tMRFZonesPhiCorr(); - forAll(*this, i) + forAll (*this, i) { operator[](i).meshPhi(MRFZonesPhiCorr); } @@ -117,7 +117,7 @@ Foam::tmp Foam::MRFZones::meshPhi() const void Foam::MRFZones::addCoriolis(fvVectorMatrix& UEqn) const { - forAll(*this, i) + forAll (*this, i) { operator[](i).addCoriolis(UEqn); } @@ -126,7 +126,7 @@ void Foam::MRFZones::addCoriolis(fvVectorMatrix& UEqn) const void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const { - forAll(*this, i) + forAll (*this, i) { operator[](i).relativeFlux(phi); } @@ -135,7 +135,7 @@ void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const void Foam::MRFZones::absoluteFlux(surfaceScalarField& phi) const { - forAll(*this, i) + forAll (*this, i) { operator[](i).absoluteFlux(phi); } @@ -148,7 +148,7 @@ void Foam::MRFZones::addCoriolis fvVectorMatrix& UEqn ) const { - forAll(*this, i) + forAll (*this, i) { operator[](i).addCoriolis(rho, UEqn); } @@ -161,7 +161,7 @@ void Foam::MRFZones::relativeFlux surfaceScalarField& phi ) const { - forAll(*this, i) + forAll (*this, i) { operator[](i).relativeFlux(rho, phi); } @@ -174,7 +174,7 @@ void Foam::MRFZones::absoluteFlux surfaceScalarField& phi ) const { - forAll(*this, i) + forAll (*this, i) { operator[](i).absoluteFlux(rho, phi); } @@ -183,7 +183,7 @@ void Foam::MRFZones::absoluteFlux void Foam::MRFZones::relativeVelocity(volVectorField& U) const { - forAll(*this, i) + forAll (*this, i) { operator[](i).relativeVelocity(U); } @@ -192,7 +192,7 @@ void Foam::MRFZones::relativeVelocity(volVectorField& U) const void Foam::MRFZones::absoluteVelocity(volVectorField& U) const { - forAll(*this, i) + forAll (*this, i) { operator[](i).absoluteVelocity(U); } @@ -201,7 +201,7 @@ void Foam::MRFZones::absoluteVelocity(volVectorField& U) const void Foam::MRFZones::correctBoundaryVelocity(volVectorField& U) const { - forAll(*this, i) + forAll (*this, i) { operator[](i).correctBoundaryVelocity(U); } @@ -237,7 +237,7 @@ Foam::tmp Foam::MRFZones::Su tmp tgradPhi = fvc::grad(phi); const volVectorField& gradPhi = tgradPhi(); - forAll(*this, i) + forAll (*this, i) { operator[](i).Su(phi, gradPhi, source); } @@ -275,7 +275,7 @@ Foam::tmp Foam::MRFZones::Su tmp tgradPhi = fvc::grad(phi); const volTensorField& gradPhi = tgradPhi(); - forAll(*this, i) + forAll (*this, i) { operator[](i).Su(phi, gradPhi, source); } diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H index 899bbf1a6..56811ceef 100644 --- a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H +++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H @@ -44,7 +44,6 @@ Description S = - (\mu \, d + \frac{\rho |U|}{2} \, f) U @f] - Since negative Darcy/Forchheimer parameters are invalid, they can be used to specify a multiplier (of the max component). @@ -106,15 +105,15 @@ class porousZone //- Coordinate system used for the zone (Cartesian) coordinateSystem coordSys_; - //- porosity of the zone (0 < porosity <= 1) + //- Porosity of the zone (0 < porosity <= 1) // Placeholder for treatment of temporal terms. // Currently unused. scalar porosity_; - //- powerLaw coefficient C0 + //- Power Law coefficient C0 scalar C0_; - //- powerLaw coefficient C1 + //- Power Law coefficient C1 scalar C1_; //- Darcy coefficient @@ -126,7 +125,7 @@ class porousZone // Private Member Functions - //- adjust negative resistance values to be multiplier of max value + //- Adjust negative resistance values to be multiplier of max value static void adjustNegativeResistance(dimensionedVector& resist); //- Power-law resistance @@ -247,7 +246,13 @@ public: return cellZoneID_; } - //- dictionary values used for the porousZone + //- Return cellZone + const cellZone& zone() const + { + return mesh_.cellZones()[cellZoneID_]; + } + + //- Dictionary values used for the porousZone const dictionary& dict() const { return dict_; diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H index 354872f93..860d1c975 100644 --- a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H +++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H @@ -92,7 +92,7 @@ class porousZones void operator=(const porousZones&); - //- modify time derivative elements + //- Modify time derivative elements template void modifyDdt(fvMatrix&) const;