From: henry <henry>
Date: Mon, 8 Jun 2009 16:41:32 +0000 (+0100)
Subject: Added support for compressible and multiphase momentum equations.
X-Git-Url: https://repo.or.cz/w/OpenFOAM-1.5.x.git/commitdiff_plain/6c35f32ce6ed2cb0819abad9aa8082037f74934e

Added support for compressible and multiphase momentum equations.
---

diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
index 8b689d0..3618d87 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,53 +32,96 @@ License
 #include "syncTools.H"
 #include "faceSet.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::MRFZone, 0);
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::MRFZone::setMRFFaces
-(
-    labelList& faceType,
-    const labelList& excludedPatchIDs
-)
+void Foam::MRFZone::setMRFFaces()
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    // Knock out coupled patches
-    forAll(patches, patchi)
+    // Type per face:
+    //  0:not in zone
+    //  1:moving with frame
+    //  2:other
+    labelList faceType(mesh_.nFaces(), 0);
+
+    // Determine faces in cell zone
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // (without constructing cells)
+
+    const labelList& own = mesh_.faceOwner();
+    const labelList& nei = mesh_.faceNeighbour();
+
+    // Cells in zone
+    boolList zoneCell(mesh_.nCells(), false);
+
+    if (cellZoneID_ != -1)
     {
-        if (patches[patchi].coupled())
+        const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
+        forAll(cellLabels, i)
         {
-            const polyPatch& pp = patches[patchi];
+            zoneCell[cellLabels[i]] = true;
+        }
+    }
 
-            forAll(pp, j)
-            {
-                label faceI = pp.start()+j;
 
-                if (faceType[faceI] == 1)
-                {
-                    faceType[faceI] = 2;
-                }
-            }
+    label nZoneFaces = 0;
+
+    for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+    {
+        if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
+        {
+            faceType[faceI] = 1;
+            nZoneFaces++;
         }
     }
 
-    // All explicitly provided exclusions
-    forAll(excludedPatchIDs, i)
+
+    labelHashSet excludedPatches(excludedPatchLabels_);
+
+    forAll(patches, patchI)
     {
-        const polyPatch& pp = patches[excludedPatchIDs[i]];
+        const polyPatch& pp = patches[patchI];
 
-        forAll(pp, j)
+        if (pp.coupled() || excludedPatches.found(patchI))
         {
-            label faceI = pp.start()+j;
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
 
-            if (faceType[faceI] == 1)
+                if (zoneCell[own[faceI]])
+                {
+                    faceType[faceI] = 2;
+                    nZoneFaces++;
+                }
+            }
+        }
+        else if (!isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
             {
-                faceType[faceI] = 2;
+                label faceI = pp.start()+i;
+
+                if (zoneCell[own[faceI]])
+                {
+                    faceType[faceI] = 1;
+                    nZoneFaces++;
+                }
             }
         }
     }
 
-    // Collect into lists per patch.
+    // Now we have for faceType:
+    //  0   : face not in cellZone
+    //  1   : internal face or normal patch face
+    //  2   : coupled patch face or excluded patch face
+
+    // Sort into lists per patch.
+
     internalFaces_.setSize(mesh_.nFaces());
     label nInternal = 0;
 
@@ -142,149 +185,43 @@ void Foam::MRFZone::setMRFFaces
         }
     }
 
-    //if (debug)
-    //{
-    //    faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
-    //    Pout<< "Writing internal faces in MRF zone to faceSet "
-    //        << internalFaces.name() << endl;
-    //    internalFaces.write();
-    //}
-    //{
-    //    faceSet MRFFaces(mesh_, "includedFaces", 100);
-    //    forAll(includedFaces_, patchi)
-    //    {
-    //        forAll(includedFaces_[patchi], i)
-    //        {
-    //            label patchFacei = includedFaces_[patchi][i];
-    //            MRFFaces.insert(patches[patchi].start()+patchFacei);
-    //        }
-    //    }
-    //    Pout<< "Writing patch faces in MRF zone to faceSet "
-    //        << MRFFaces.name() << endl;
-    //    MRFFaces.write();
-    //}
-    //{
-    //    faceSet excludedFaces(mesh_, "excludedFaces", 100);
-    //    forAll(excludedFaces_, patchi)
-    //    {
-    //        forAll(excludedFaces_[patchi], i)
-    //        {
-    //            label patchFacei = excludedFaces_[patchi][i];
-    //            excludedFaces.insert(patches[patchi].start()+patchFacei);
-    //        }
-    //    }
-    //    Pout<< "Writing faces in MRF zone with special handling to faceSet "
-    //        << excludedFaces.name() << endl;
-    //    excludedFaces.write();
-    //}
-}
-
-
-void Foam::MRFZone::setMRFFaces()
-{
-    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-    // Type per face:
-    //  0:not in zone
-    //  1:moving with frame
-    //  2:other
-    labelList faceType(mesh_.nFaces(), 0);
-
-    bool faceZoneFound = (faceZoneID_ != -1);
-    reduce(faceZoneFound, orOp<bool>());
 
-    if (faceZoneFound)
+    if (debug)
     {
-        // Explicitly provided faces.
-        if (faceZoneID_ != -1)
+        faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
+        Pout<< "Writing " << internalFaces.size()
+            << " internal faces in MRF zone to faceSet "
+            << internalFaces.name() << endl;
+        internalFaces.write();
+
+        faceSet MRFFaces(mesh_, "includedFaces", 100);
+        forAll(includedFaces_, patchi)
         {
-            const labelList& zoneFaces = mesh_.faceZones()[faceZoneID_];
-
-            forAll(zoneFaces, i)
-            {
-                faceType[zoneFaces[i]] = 1;
-            }
-
-            if (allPatchesMove_)
+            forAll(includedFaces_[patchi], i)
             {
-                // Explicitly provided patches that do not move
-                setMRFFaces(faceType, patchLabels_);
-            }
-            else
-            {
-                setMRFFaces(faceType, labelList(0));
+                label patchFacei = includedFaces_[patchi][i];
+                MRFFaces.insert(patches[patchi].start()+patchFacei);
             }
         }
-    }
-    else
-    {
-        // Determine faces in cell zone
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        // (without constructing cells)
-
-        const labelList& own = mesh_.faceOwner();
-        const labelList& nei = mesh_.faceNeighbour();
-
-        // Cells in zone
-        boolList zoneCell(mesh_.nCells(), false);
+        Pout<< "Writing " << MRFFaces.size()
+            << " patch faces in MRF zone to faceSet "
+            << MRFFaces.name() << endl;
+        MRFFaces.write();
 
-        if (cellZoneID_ != -1)
+        faceSet excludedFaces(mesh_, "excludedFaces", 100);
+        forAll(excludedFaces_, patchi)
         {
-            const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
-            forAll(cellLabels, i)
+            forAll(excludedFaces_[patchi], i)
             {
-                zoneCell[cellLabels[i]] = true;
+                label patchFacei = excludedFaces_[patchi][i];
+                excludedFaces.insert(patches[patchi].start()+patchFacei);
             }
         }
-
-
-        label nZoneFaces = 0;
-
-        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
-        {
-            if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
-            {
-                faceType[faceI] = 1;
-                nZoneFaces++;
-            }
-        }
-        forAll(patches, patchI)
-        {
-            const polyPatch& pp = patches[patchI];
-
-            if (!isA<emptyPolyPatch>(pp))
-            {
-                forAll(pp, i)
-                {
-                    label faceI = pp.start()+i;
-
-                    if (zoneCell[own[faceI]])
-                    {
-                        faceType[faceI] = 1;
-                        nZoneFaces++;
-                    }
-                }
-            }
-        }
-        syncTools::syncFaceList(mesh_, faceType, maxEqOp<label>(), false);
-
-
-        Info<< nl
-            << "MRFZone " << name_ << " : did not find a faceZone; using "
-            << returnReduce(nZoneFaces, sumOp<label>())
-            << " faces from the cellZone instead." << endl;
-
-
-        if (allPatchesMove_)
-        {
-            // Explicitly provided excluded patches
-            setMRFFaces(faceType, patchLabels_);
-        }
-        else
-        {
-            setMRFFaces(faceType, labelList(0));
-        }
-    }    
+        Pout<< "Writing " << excludedFaces.size()
+            << " faces in MRF zone with special handling to faceSet "
+            << excludedFaces.name() << endl;
+        excludedFaces.write();
+    }
 }
 
 
@@ -296,36 +233,41 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
     name_(is),
     dict_(is),
     cellZoneID_(mesh_.cellZones().findZoneID(name_)),
-    faceZoneID_(mesh_.faceZones().findZoneID(name_)),
-    allPatchesMove_(dict_.found("nonRotatingPatches")),
-    patchNames_
+    excludedPatchNames_
     (
-        allPatchesMove_
-      ? dict_.lookup("nonRotatingPatches")
-      : dict_.lookup("patches")
+        dict_.lookupOrDefault("nonRotatingPatches", wordList(0))
     ),
     origin_(dict_.lookup("origin")),
     axis_(dict_.lookup("axis")),
     omega_(dict_.lookup("omega")),
     Omega_("Omega", omega_*axis_)
 {
+    if (dict_.found("patches"))
+    {
+        WarningIn("MRFZone(const fvMesh&, Istream&)")
+            << "Ignoring entry 'patches'\n"
+            << "    By default all patches within the rotating region rotate.\n"
+            << "    Optionally supply excluded patches using 'nonRotatingPatches'."
+            << endl;
+    }
+
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
     axis_ = axis_/mag(axis_);
     Omega_ = omega_*axis_;
 
-    patchLabels_.setSize(patchNames_.size());
+    excludedPatchLabels_.setSize(excludedPatchNames_.size());
 
-    forAll(patchNames_, i)
+    forAll(excludedPatchNames_, i)
     {
-        patchLabels_[i] = patches.findPatchID(patchNames_[i]);
+        excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]);
 
-        if (patchLabels_[i] == -1)
+        if (excludedPatchLabels_[i] == -1)
         {
             FatalErrorIn
             (
                 "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
-            )   << "cannot find MRF patch " << patchNames_[i]
+            )   << "cannot find MRF patch " << excludedPatchNames_[i]
                 << exit(FatalError);
         }
     }
@@ -364,7 +306,71 @@ void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
 
     forAll(cells, i)
     {
-        Usource[cells[i]] -= V[cells[i]]*(Omega ^ U[cells[i]]);
+        label celli = cells[i];
+        Usource[celli] -= V[celli]*(Omega ^ U[celli]);
+    }
+}
+
+
+void Foam::MRFZone::addCoriolis
+(
+    const volScalarField& rho,
+    fvVectorMatrix& UEqn
+) const
+{
+    if (cellZoneID_ == -1)
+    {
+        return;
+    }
+
+    const labelList& cells = mesh_.cellZones()[cellZoneID_];
+    const scalarField& V = mesh_.V();
+    vectorField& Usource = UEqn.source();
+    const vectorField& U = UEqn.psi();
+    const vector& Omega = Omega_.value();
+
+    forAll(cells, i)
+    {
+        label celli = cells[i];
+        Usource[celli] -= V[celli]*rho[celli]*(Omega ^ U[celli]);
+    }
+}
+
+
+void Foam::MRFZone::relativeVelocity(volVectorField& U) const
+{
+    const volVectorField& C = mesh_.C();
+
+    const vector& origin = origin_.value();
+    const vector& Omega = Omega_.value();
+
+    const labelList& cells = mesh_.cellZones()[cellZoneID_];
+
+    forAll(cells, i)
+    {
+        label celli = cells[i];
+        U[celli] -= (Omega ^ (C[celli] - origin));
+    }
+
+    // Included patches
+    forAll(includedFaces_, patchi)
+    {
+        forAll(includedFaces_[patchi], i)
+        {
+            label patchFacei = includedFaces_[patchi][i];
+            U.boundaryField()[patchi][patchFacei] = vector::zero;
+        }
+    }
+
+    // Excluded patches
+    forAll(excludedFaces_, patchi)
+    {
+        forAll(excludedFaces_[patchi], i)
+        {
+            label patchFacei = excludedFaces_[patchi][i];
+            U.boundaryField()[patchi][patchFacei] -=
+                (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
+        }
     }
 }
 
@@ -410,6 +416,49 @@ void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
 }
 
 
+void Foam::MRFZone::absoluteFlux(surfaceScalarField& phi) const
+{
+    const surfaceVectorField& Cf = mesh_.Cf();
+    const surfaceVectorField& Sf = mesh_.Sf();
+
+    const vector& origin = origin_.value();
+    const vector& Omega = Omega_.value();
+
+    // Internal faces
+    forAll(internalFaces_, i)
+    {
+        label facei = internalFaces_[i];
+        phi[facei] += (Omega ^ (Cf[facei] - origin)) & Sf[facei];
+    }
+
+    // Included patches
+    forAll(includedFaces_, patchi)
+    {
+        forAll(includedFaces_[patchi], i)
+        {
+            label patchFacei = includedFaces_[patchi][i];
+
+            phi.boundaryField()[patchi][patchFacei] +=
+                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+              & Sf.boundaryField()[patchi][patchFacei];
+        }
+    }
+
+    // Excluded patches
+    forAll(excludedFaces_, patchi)
+    {
+        forAll(excludedFaces_[patchi], i)
+        {
+            label patchFacei = excludedFaces_[patchi][i];
+
+            phi.boundaryField()[patchi][patchFacei] +=
+                (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
+              & Sf.boundaryField()[patchi][patchFacei];
+        }
+    }
+}
+
+
 void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
 {
     const vector& origin = origin_.value();
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
index 5926e87..d5554a5 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,7 +26,7 @@ Class
     Foam::MRFZone
 
 Description
-    MRF zone definition based on cell zone and optional face zone and parameters
+    MRF zone definition based on cell zone and parameters
     obtained from a control dictionary constructed from the given stream.
 
     The rotation of the MRF region is defined by an origin and axis of
@@ -74,15 +74,8 @@ class MRFZone
 
         label cellZoneID_;
 
-        //- label of face zone with faces on outside of cell zone.
-        label faceZoneID_;
-
-        //- Do patches move with frame (true) or are explicitly provided (false,
-        //  old behaviour)
-        Switch allPatchesMove_;
-
-        const wordList patchNames_;
-        labelList patchLabels_;
+        const wordList excludedPatchNames_;
+        labelList excludedPatchLabels_;
 
         //- Internal faces that are part of MRF
         labelList internalFaces_;
@@ -102,16 +95,8 @@ class MRFZone
     // Private Member Functions
 
         //- Divide faces in frame according to patch
-        void setMRFFaces
-        (
-            labelList& faceType,
-            const labelList& excludedPatchIDs
-        );
-
-        //- Divide faces in frame according to patch
         void setMRFFaces();
 
-
         //- Disallow default bitwise copy construct
         MRFZone(const MRFZone&);
 
@@ -121,6 +106,10 @@ class MRFZone
 
 public:
 
+    // Declare name of the class and its debug switch
+    ClassName("MRFZone");
+
+
     // Constructors
 
         //- Construct from fvMesh and Istream
@@ -165,9 +154,18 @@ public:
         //- Add the Coriolis force contribution to the momentum equation
         void addCoriolis(fvVectorMatrix& UEqn) const;
 
+        //- Add the Coriolis force contribution to the momentum equation
+        void addCoriolis(const volScalarField& rho, fvVectorMatrix& UEqn) const;
+
+        //- Make the given absolute velocity relative within the MRF region
+        void relativeVelocity(volVectorField& U) const;
+
         //- Make the given absolute flux relative within the MRF region
         void relativeFlux(surfaceScalarField& phi) const;
 
+        //- Make the given relative flux absolute within the MRF region
+        void absoluteFlux(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/MRFZones.C b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C
index 4d7389b..5ef8b3e 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.C
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,6 +65,28 @@ void Foam::MRFZones::addCoriolis(fvVectorMatrix& UEqn) const
 }
 
 
+void Foam::MRFZones::addCoriolis
+(
+    const volScalarField& rho,
+    fvVectorMatrix& UEqn
+) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).addCoriolis(rho, UEqn);
+    }
+}
+
+
+void Foam::MRFZones::relativeVelocity(volVectorField& U) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).relativeVelocity(U);
+    }
+}
+
+
 void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const
 {
     forAll(*this, i)
@@ -74,6 +96,15 @@ void Foam::MRFZones::relativeFlux(surfaceScalarField& phi) const
 }
 
 
+void Foam::MRFZones::absoluteFlux(surfaceScalarField& phi) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).absoluteFlux(phi);
+    }
+}
+
+
 void Foam::MRFZones::correctBoundaryVelocity(volVectorField& U) const
 {
     forAll(*this, i)
diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZones.H b/src/finiteVolume/cfdTools/general/MRF/MRFZones.H
index 5965b86..6e97990 100644
--- a/src/finiteVolume/cfdTools/general/MRF/MRFZones.H
+++ b/src/finiteVolume/cfdTools/general/MRF/MRFZones.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,9 +76,18 @@ public:
         //- Add the Coriolis force contribution to the momentum equation
         void addCoriolis(fvVectorMatrix& UEqn) const;
 
+        //- Add the Coriolis force contribution to the momentum equation
+        void addCoriolis(const volScalarField& rho, fvVectorMatrix& UEqn) const;
+
+        //- Make the given absolute velocity relative within the MRF region
+        void relativeVelocity(volVectorField& U) const;
+
         //- Make the given absolute flux relative within the MRF region
         void relativeFlux(surfaceScalarField& phi) const;
 
+        //- Make the given relative flux absolute within the MRF region
+        void absoluteFlux(surfaceScalarField& phi) const;
+
         //- Correct the boundary velocity for the roation of the MRF region
         void correctBoundaryVelocity(volVectorField& U) const;
 };