Updates to MRF handling for block-coupled solver

applications/solvers/coupled/MRFPorousFoam/UEqn.H

@@ -9,13 +9,10 @@
       + turbulence->divDevReff()
     );
 
-    // Add MRF sources
+    // Add MRF and porous sources implicitly.  HJ, 18/Nov/2017
-    mrfZones.addCoriolis(UEqn);
-
-    // Add porous sources
     tmp<volTensorField> tTU;
 
-    if (addPorosity)
+    if (addMRF || addPorosity)
     {
         tTU = tmp<volTensorField>
         (
@@ -35,6 +32,10 @@
         );
         volTensorField& TU = tTU();
 
+        // Add implicit MRF source as a Hodge dual of the rotational velocity
+        TU += *mrfZones.omega();
+        
+        // Add implicit resistance
         pZones.addResistance(UEqn, TU);
 
         trTU = inv(TU + tensor(I)*UEqn.A());
@@ -52,7 +53,7 @@
     // Insert momentum equation
     UpEqn.insertEquation(0, UEqn);
 
-    if (addPorosity)
+    if (addMRF || addPorosity)
     {
         // Manually over-ride the 3x3 block to handle the off-diagonal
         // part of the Ap coefficient
@@ -62,11 +63,43 @@
 
         const scalarField& V = mesh.V().field();
 
-        // Note: insertion should only happen in porous cell zones
+        // Note: insertion should only happen in MRF and porous cell zones
         // HJ, 14/Mar/2016
 
+        // Warning.  Possible problem with a zone which is both MRF and porous
+        // The solution is to do the loop below everywhere
+        // Currently only inserting zones for efficiency.  HJ, 18/Nov/2017
         register label cellI;
 
+        forAll (mrfZones, mrfZoneI)
+        {
+            const labelList& curZoneCells = mrfZones[mrfZoneI].zone();
+
+            // Loop over all cells in the zone
+            forAll (curZoneCells, zcI)
+            {
+                cellI = curZoneCells[zcI];
+
+                const scalar& cellV = V[cellI];
+
+                const tensor& cellTU = TUIn[cellI];
+
+                CoeffField<vector4>::squareType& cellDD = DD[cellI];
+
+                cellDD(0, 0) += cellV*cellTU.xx();
+                cellDD(0, 1) += cellV*cellTU.xy();
+                cellDD(0, 2) += cellV*cellTU.xz();
+
+                cellDD(1, 0) += cellV*cellTU.yx();
+                cellDD(1, 1) += cellV*cellTU.yy();
+                cellDD(2, 2) += cellV*cellTU.yz();
+
+                cellDD(2, 0) += cellV*cellTU.zx();
+                cellDD(2, 1) += cellV*cellTU.zy();
+                cellDD(2, 2) += cellV*cellTU.zz();
+            }
+        }
+
         forAll (pZones, pZoneI)
         {
             const labelList& curZoneCells = pZones[pZoneI].zone();

applications/solvers/coupled/MRFPorousFoam/createMRF.H

@@ -1,2 +1,9 @@
     MRFZones mrfZones(mesh);
     mrfZones.correctBoundaryVelocity(U);
+
+    bool addMRF = false;
+
+    if (!mrfZones.empty())
+    {
+        addMRF = true;
+    }

applications/solvers/coupled/MRFPorousFoam/pEqn.H

@@ -2,7 +2,7 @@
 tmp<fvScalarMatrix> tpEqn;
 tmp<surfaceScalarField> tpresSource;
 
-if (addPorosity)
+if (addMRF || addPorosity)
 {
     // Collect pressure source with tensorial 1/Ap
     const volTensorField& rTU = trTU();

src/finiteVolume/cfdTools/general/MRF/MRFZone.C

@@ -581,6 +581,20 @@ void Foam::MRFZone::addCoriolis
 }
 
 
+void Foam::MRFZone::addOmega(volVectorField& omg) const
+{
+    const vector rotVel = Omega();
+
+    // Set omega in all cells of the rotating cell zone
+    const labelList& cells = mesh_.cellZones()[cellZoneID_];
+
+    forAll (cells, i)
+    {
+        omg[cells[i]] = rotVel;
+    }
+}
+
+
 void Foam::MRFZone::relativeVelocity(volVectorField& U) const
 {
     const volVectorField& C = mesh_.C();

src/finiteVolume/cfdTools/general/MRF/MRFZone.H

@@ -111,7 +111,7 @@ class MRFZone
 
         //- Time for which mesh velocity is calculated
         mutable scalar meshVelTime_;
-    
+
 
     // Private Static Data
 
@@ -203,6 +203,18 @@ public:
                 return name_;
             }
 
+            //- Return cellZone number
+            label zoneId() const
+            {
+                return cellZoneID_;
+            }
+
+            //- Return cellZone
+            const cellZone& zone() const
+            {
+                return mesh_.cellZones()[cellZoneID_];
+            }
+
             //- Return the MRFZone dictionary
             const dictionary& dict() const
             {
@@ -235,7 +247,7 @@ public:
         {
             // Clear mesh velocity
             deleteDemandDrivenData(meshVelocityPtr_);
-            
+
             // Only updates face addressing
             setMRFFaces();
         }
@@ -250,6 +262,9 @@ public:
             fvVectorMatrix& UEqn
         ) const;
 
+        //- Add omega within the MRF region
+        void addOmega(volVectorField& omg) const;
+
         //- Make the given absolute velocity relative within the MRF region
         void relativeVelocity(volVectorField& U) const;
 
@@ -276,7 +291,7 @@ public:
             surfaceScalarField& phi
         ) const;
 
-        //- Compute the pseudo mesh phi of the MRF region
+        //- Compute pseudo mesh phi for the MRF region
         void meshPhi(surfaceScalarField& phi) const;
 
         //- Correct the boundary velocity for the roation of the MRF region

src/finiteVolume/cfdTools/general/MRF/MRFZones.C

@@ -57,6 +57,35 @@ Foam::MRFZones::MRFZones(const fvMesh& mesh)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::tmp<Foam::volVectorField> Foam::MRFZones::omega() const
+{
+    tmp<volVectorField> tMRFZonesOmega
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                "MRFZonesOmega",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector("zero", dimless/dimTime, vector::zero)
+        )
+    );
+    volVectorField& MRFZonesOmega = tMRFZonesOmega();
+
+    forAll (*this, i)
+    {
+        operator[](i).addOmega(MRFZonesOmega);
+    }
+
+    return tMRFZonesOmega;
+}
+
+
 Foam::tmp<Foam::surfaceScalarField> Foam::MRFZones::fluxCorrection() const
 {
     tmp<surfaceScalarField> tMRFZonesPhiCorr
@@ -88,7 +117,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::MRFZones::fluxCorrection() const
 
 Foam::tmp<Foam::surfaceScalarField> Foam::MRFZones::meshPhi() const
 {
-    tmp<surfaceScalarField> tMRFZonesPhiCorr
+    tmp<surfaceScalarField> tMRFZonesFaceU
     (
         new surfaceScalarField
         (
@@ -104,14 +133,14 @@ Foam::tmp<Foam::surfaceScalarField> Foam::MRFZones::meshPhi() const
             dimensionedScalar("zero", dimVolume/dimTime, 0)
         )
     );
-    surfaceScalarField& MRFZonesPhiCorr = tMRFZonesPhiCorr();
+    surfaceScalarField& MRFZonesFaceU = tMRFZonesFaceU();
 
     forAll (*this, i)
     {
-        operator[](i).meshPhi(MRFZonesPhiCorr);
+        operator[](i).meshPhi(MRFZonesFaceU);
     }
 
-    return tMRFZonesPhiCorr;
+    return tMRFZonesFaceU;
 }
 
 

src/finiteVolume/cfdTools/general/MRF/MRFZones.H

@@ -77,12 +77,16 @@ public:
 
     // Member Functions
 
+        //- Return cell centre omega over all zones
+        tmp<volVectorField> omega() const;
+
         //- Return raw correction flux
         tmp<surfaceScalarField> fluxCorrection() const;
 
         //- Return pseudo mesh flux
         tmp<surfaceScalarField> meshPhi() const;
 
+
         // Incompressible MRF
 
             //- Add the Coriolis force contribution to the momentum equation