From 78a9142214076cf5533254883d777328cbc97ab1 Mon Sep 17 00:00:00 2001 From: Hrvoje Jasak Date: Tue, 24 Jan 2017 16:44:22 +0000 Subject: [PATCH] Bugfix: MRF zone flux check --- .../cfdTools/general/MRF/MRFZone.C | 60 +++++++++------- .../cfdTools/general/MRF/MRFZone.H | 2 +- .../cfdTools/general/MRF/MRFZoneTemplates.C | 70 +++++++++++-------- .../cfdTools/general/MRF/MRFZones.C | 10 +-- .../cfdTools/general/MRF/MRFZones.H | 4 +- 5 files changed, 86 insertions(+), 60 deletions(-) diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C index 0390307f2..9633fff7a 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C @@ -139,13 +139,15 @@ void Foam::MRFZone::setMRFFaces() labelList nIncludedFaces(patches.size(), 0); labelList nExcludedFaces(patches.size(), 0); + register label faceI; + forAll (patches, patchi) { const polyPatch& pp = patches[patchi]; forAll (pp, patchFaceI) { - label faceI = pp.start() + patchFaceI; + faceI = pp.start() + patchFaceI; if (faceType[faceI] == 1) { @@ -174,7 +176,7 @@ void Foam::MRFZone::setMRFFaces() forAll (pp, patchFaceI) { - label faceI = pp.start() + patchFaceI; + faceI = pp.start() + patchFaceI; if (faceType[faceI] == 1) { @@ -187,7 +189,6 @@ void Foam::MRFZone::setMRFFaces() } } - if (debug) { faceSet internalFaces @@ -489,42 +490,51 @@ void Foam::MRFZone::absoluteFlux void Foam::MRFZone::meshPhi ( - surfaceVectorField& phi + surfaceScalarField& phi ) const { const surfaceVectorField& Cf = mesh_.Cf(); + const surfaceVectorField& Sf = mesh_.Sf(); const vector& origin = origin_.value(); const vector rotVel = Omega(); // Internal faces + const vectorField& CfIn = Cf.internalField(); + const vectorField& SfIn = Sf.internalField(); + + register label faceI, patchFaceI; + forAll (internalFaces_, i) { - label facei = internalFaces_[i]; - phi[facei] = (rotVel ^ (Cf[facei] - origin)); + faceI = internalFaces_[i]; + phi[faceI] = SfIn[faceI] & (rotVel ^ (CfIn[faceI] - origin)); } // Included patches - forAll (includedFaces_, patchi) - { - forAll (includedFaces_[patchi], i) - { - label patchFaceI = includedFaces_[patchi][i]; - phi.boundaryField()[patchi][patchFaceI] = - (rotVel ^ (Cf.boundaryField()[patchi][patchFaceI] - origin)); + forAll (includedFaces_, patchI) + { + forAll (includedFaces_[patchI], i) + { + patchFaceI = includedFaces_[patchI][i]; + + phi.boundaryField()[patchI][patchFaceI] = + Sf.boundaryField()[patchI][patchFaceI] + & (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin)); } } // Excluded patches - forAll (excludedFaces_, patchi) + forAll (excludedFaces_, patchI) { - forAll (excludedFaces_[patchi], i) + forAll (excludedFaces_[patchI], i) { - label patchFaceI = excludedFaces_[patchi][i]; + patchFaceI = excludedFaces_[patchI][i]; - phi.boundaryField()[patchi][patchFaceI] = - (rotVel ^ (Cf.boundaryField()[patchi][patchFaceI] - origin)); + phi.boundaryField()[patchI][patchFaceI] = + Sf.boundaryField()[patchI][patchFaceI] + & (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin)); } } } @@ -534,21 +544,23 @@ void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const const vector& origin = origin_.value(); const vector rotVel = Omega(); + register label patchFaceI; + // Included patches - forAll (includedFaces_, patchi) + forAll (includedFaces_, patchI) { - const vectorField& patchC = mesh_.Cf().boundaryField()[patchi]; + const vectorField& patchC = mesh_.Cf().boundaryField()[patchI]; - vectorField pfld(U.boundaryField()[patchi]); + vectorField pfld(U.boundaryField()[patchI]); - forAll (includedFaces_[patchi], i) + forAll (includedFaces_[patchI], i) { - label patchFaceI = includedFaces_[patchi][i]; + patchFaceI = includedFaces_[patchI][i]; pfld[patchFaceI] = (rotVel ^ (patchC[patchFaceI] - origin)); } - U.boundaryField()[patchi] == pfld; + U.boundaryField()[patchI] == pfld; } } diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H index 954bc85b1..0121eed2a 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H @@ -221,7 +221,7 @@ public: ) const; //- Compute the pseudo mesh phi of the MRF region - void meshPhi(surfaceVectorField& phi) const; + void meshPhi(surfaceScalarField& phi) const; //- Correct the boundary velocity for the roation of the MRF region void correctBoundaryVelocity(volVectorField& U) const; diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C index 0ae07364a..a6eabcbfa 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C @@ -44,34 +44,42 @@ void Foam::MRFZone::relativeRhoFlux const vector rotVel = Omega(); // Internal faces + const vectorField& CfIn = Cf.internalField(); + const vectorField& SfIn = Sf.internalField(); + + register label faceI, patchFaceI; + forAll (internalFaces_, i) { - label facei = internalFaces_[i]; - phi[facei] -= rho[facei]*(rotVel ^ (Cf[facei] - origin)) & Sf[facei]; + faceI = internalFaces_[i]; + + phi[faceI] -= + rho[faceI]*(rotVel ^ (CfIn[faceI] - origin)) & SfIn[faceI]; } - // Included patches - forAll (includedFaces_, patchi) + // Included patches: reset the flux to exactly zero to avoid + // round-off issues + forAll (includedFaces_, patchI) { - forAll (includedFaces_[patchi], i) + forAll (includedFaces_[patchI], i) { - label patchFacei = includedFaces_[patchi][i]; + patchFaceI = includedFaces_[patchI][i]; - phi.boundaryField()[patchi][patchFacei] = 0.0; + phi.boundaryField()[patchI][patchFaceI] = 0.0; } } // Excluded patches - forAll (excludedFaces_, patchi) + forAll (excludedFaces_, patchI) { - forAll (excludedFaces_[patchi], i) + forAll (excludedFaces_[patchI], i) { - label patchFacei = excludedFaces_[patchi][i]; + patchFaceI = excludedFaces_[patchI][i]; - phi.boundaryField()[patchi][patchFacei] -= - rho.boundaryField()[patchi][patchFacei] - *(rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)) - & Sf.boundaryField()[patchi][patchFacei]; + phi.boundaryField()[patchI][patchFaceI] -= + rho.boundaryField()[patchI][patchFaceI] + *(rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin)) + & Sf.boundaryField()[patchI][patchFaceI]; } } } @@ -91,35 +99,41 @@ void Foam::MRFZone::absoluteRhoFlux const vector rotVel = Omega(); // Internal faces + const vectorField& CfIn = Cf.internalField(); + const vectorField& SfIn = Sf.internalField(); + + register label faceI, patchFaceI; + forAll (internalFaces_, i) { - label facei = internalFaces_[i]; - phi[facei] += (rotVel ^ (Cf[facei] - origin)) & Sf[facei]; + faceI = internalFaces_[i]; + + phi[faceI] += (rotVel ^ (CfIn[faceI] - origin)) & SfIn[faceI]; } // Included patches - forAll (includedFaces_, patchi) + forAll (includedFaces_, patchI) { - forAll (includedFaces_[patchi], i) + forAll (includedFaces_[patchI], i) { - label patchFacei = includedFaces_[patchi][i]; + patchFaceI = includedFaces_[patchI][i]; - phi.boundaryField()[patchi][patchFacei] += - (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)) - & Sf.boundaryField()[patchi][patchFacei]; + phi.boundaryField()[patchI][patchFaceI] += + (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin)) + & Sf.boundaryField()[patchI][patchFaceI]; } } // Excluded patches - forAll (excludedFaces_, patchi) + forAll (excludedFaces_, patchI) { - forAll (excludedFaces_[patchi], i) + forAll (excludedFaces_[patchI], i) { - label patchFacei = excludedFaces_[patchi][i]; + patchFaceI = excludedFaces_[patchI][i]; - phi.boundaryField()[patchi][patchFacei] += - (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)) - & Sf.boundaryField()[patchi][patchFacei]; + phi.boundaryField()[patchI][patchFaceI] += + (rotVel ^ (Cf.boundaryField()[patchI][patchFaceI] - origin)) + & Sf.boundaryField()[patchI][patchFaceI]; } } } diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C index e5f29c690..f9cb77151 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C @@ -86,11 +86,11 @@ Foam::tmp Foam::MRFZones::fluxCorrection() const } -Foam::tmp Foam::MRFZones::meshPhi() const +Foam::tmp Foam::MRFZones::meshPhi() const { - tmp tMRFZonesPhiCorr + tmp tMRFZonesPhiCorr ( - new surfaceVectorField + new surfaceScalarField ( IOobject ( @@ -101,10 +101,10 @@ Foam::tmp Foam::MRFZones::meshPhi() const IOobject::NO_WRITE ), mesh_, - dimensionedVector("zero", dimVelocity, vector::zero) + dimensionedScalar("zero", dimVolume/dimTime, 0) ) ); - surfaceVectorField& MRFZonesPhiCorr = tMRFZonesPhiCorr(); + surfaceScalarField& MRFZonesPhiCorr = tMRFZonesPhiCorr(); forAll (*this, i) { diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZones.H b/src/finiteVolume/cfdTools/general/MRF/MRFZones.H index ff683525d..082fc61f2 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.H +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.H @@ -81,7 +81,7 @@ public: tmp fluxCorrection() const; //- Return pseudo mesh flux - tmp meshPhi() const; + tmp meshPhi() const; // Incompressible MRF @@ -134,7 +134,7 @@ public: void absoluteVelocity(volVectorField& Urel) const; - //- Correct the boundary velocity for the roation of the MRF region + //- Correct the boundary velocity for the rotation of the MRF region void correctBoundaryVelocity(volVectorField& U) const; //- Compute source term for volScalarFields