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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "MRFZone.H"
  28 #include "fvMesh.H"
  29 #include "volFields.H"
  30 #include "surfaceFields.H"
  31 #include "fvMatrices.H"
  32 #include "syncTools.H"
  33 #include "faceSet.H"
  34 #include "geometricOneField.H"
  35
  36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  37
  38 defineTypeNameAndDebug(Foam::MRFZone, 0);
  39
  40
  41 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  42
  43 void Foam::MRFZone::setMRFFaces()
  44 {
  45     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
  46
  47     // Type per face:
  48     //  0:not in zone
  49     //  1:moving with frame
  50     //  2:other
  51     labelList faceType(mesh_.nFaces(), 0);
  52
  53     // Determine faces in cell zone
  54     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  55     // (without constructing cells)
  56
  57     const labelList& own = mesh_.faceOwner();
  58     const labelList& nei = mesh_.faceNeighbour();
  59
  60     // Cells in zone
  61     boolList zoneCell(mesh_.nCells(), false);
  62
  63     if (cellZoneID_ != -1)
  64     {
  65         const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
  66         forAll(cellLabels, i)
  67         {
  68             zoneCell[cellLabels[i]] = true;
  69         }
  70     }
  71
  72
  73     label nZoneFaces = 0;
  74
  75     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
  76     {
  77         if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
  78         {
  79             faceType[faceI] = 1;
  80             nZoneFaces++;
  81         }
  82     }
  83
  84
  85     labelHashSet excludedPatches(excludedPatchLabels_);
  86
  87     forAll(patches, patchI)
  88     {
  89         const polyPatch& pp = patches[patchI];
  90
  91         if (pp.coupled() || excludedPatches.found(patchI))
  92         {
  93             forAll(pp, i)
  94             {
  95                 label faceI = pp.start()+i;
  96
  97                 if (zoneCell[own[faceI]])
  98                 {
  99                     faceType[faceI] = 2;
 100                     nZoneFaces++;
 101                 }
 102             }
 103         }
 104         else if (!isA<emptyPolyPatch>(pp))
 105         {
 106             forAll(pp, i)
 107             {
 108                 label faceI = pp.start()+i;
 109
 110                 if (zoneCell[own[faceI]])
 111                 {
 112                     faceType[faceI] = 1;
 113                     nZoneFaces++;
 114                 }
 115             }
 116         }
 117     }
 118
 119     // Now we have for faceType:
 120     //  0   : face not in cellZone
 121     //  1   : internal face or normal patch face
 122     //  2   : coupled patch face or excluded patch face
 123
 124     // Sort into lists per patch.
 125
 126     internalFaces_.setSize(mesh_.nFaces());
 127     label nInternal = 0;
 128
 129     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 130     {
 131         if (faceType[faceI] == 1)
 132         {
 133             internalFaces_[nInternal++] = faceI;
 134         }
 135     }
 136     internalFaces_.setSize(nInternal);
 137
 138     labelList nIncludedFaces(patches.size(), 0);
 139     labelList nExcludedFaces(patches.size(), 0);
 140
 141     forAll(patches, patchi)
 142     {
 143         const polyPatch& pp = patches[patchi];
 144
 145         forAll(pp, patchFacei)
 146         {
 147             label faceI = pp.start()+patchFacei;
 148
 149             if (faceType[faceI] == 1)
 150             {
 151                 nIncludedFaces[patchi]++;
 152             }
 153             else if (faceType[faceI] == 2)
 154             {
 155                 nExcludedFaces[patchi]++;
 156             }
 157         }
 158     }
 159
 160     includedFaces_.setSize(patches.size());
 161     excludedFaces_.setSize(patches.size());
 162     forAll(nIncludedFaces, patchi)
 163     {
 164         includedFaces_[patchi].setSize(nIncludedFaces[patchi]);
 165         excludedFaces_[patchi].setSize(nExcludedFaces[patchi]);
 166     }
 167     nIncludedFaces = 0;
 168     nExcludedFaces = 0;
 169
 170     forAll(patches, patchi)
 171     {
 172         const polyPatch& pp = patches[patchi];
 173
 174         forAll(pp, patchFacei)
 175         {
 176             label faceI = pp.start()+patchFacei;
 177
 178             if (faceType[faceI] == 1)
 179             {
 180                 includedFaces_[patchi][nIncludedFaces[patchi]++] = patchFacei;
 181             }
 182             else if (faceType[faceI] == 2)
 183             {
 184                 excludedFaces_[patchi][nExcludedFaces[patchi]++] = patchFacei;
 185             }
 186         }
 187     }
 188
 189
 190     if (debug)
 191     {
 192         faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
 193         Pout<< "Writing " << internalFaces.size()
 194             << " internal faces in MRF zone to faceSet "
 195             << internalFaces.name() << endl;
 196         internalFaces.write();
 197
 198         faceSet MRFFaces(mesh_, "includedFaces", 100);
 199         forAll(includedFaces_, patchi)
 200         {
 201             forAll(includedFaces_[patchi], i)
 202             {
 203                 label patchFacei = includedFaces_[patchi][i];
 204                 MRFFaces.insert(patches[patchi].start()+patchFacei);
 205             }
 206         }
 207         Pout<< "Writing " << MRFFaces.size()
 208             << " patch faces in MRF zone to faceSet "
 209             << MRFFaces.name() << endl;
 210         MRFFaces.write();
 211
 212         faceSet excludedFaces(mesh_, "excludedFaces", 100);
 213         forAll(excludedFaces_, patchi)
 214         {
 215             forAll(excludedFaces_[patchi], i)
 216             {
 217                 label patchFacei = excludedFaces_[patchi][i];
 218                 excludedFaces.insert(patches[patchi].start()+patchFacei);
 219             }
 220         }
 221         Pout<< "Writing " << excludedFaces.size()
 222             << " faces in MRF zone with special handling to faceSet "
 223             << excludedFaces.name() << endl;
 224         excludedFaces.write();
 225     }
 226 }
 227
 228
 229 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 230
 231 Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
 232 :
 233     mesh_(mesh),
 234     name_(is),
 235     dict_(is),
 236     cellZoneID_(mesh_.cellZones().findZoneID(name_)),
 237     excludedPatchNames_
 238     (
 239         dict_.lookupOrDefault("nonRotatingPatches", wordList(0))
 240     ),
 241     origin_(dict_.lookup("origin")),
 242     axis_(dict_.lookup("axis")),
 243     omega_(dict_.lookup("omega")),
 244     Omega_("Omega", omega_*axis_)
 245 {
 246     if (dict_.found("patches"))
 247     {
 248         WarningIn("MRFZone(const fvMesh&, Istream&)")
 249             << "Ignoring entry 'patches'\n"
 250             << "    By default all patches within the rotating region rotate.\n"
 251             << "    Optionally supply excluded patches using 'nonRotatingPatches'."
 252             << endl;
 253     }
 254
 255     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 256
 257     axis_ = axis_/mag(axis_);
 258     Omega_ = omega_*axis_;
 259
 260     excludedPatchLabels_.setSize(excludedPatchNames_.size());
 261
 262     forAll(excludedPatchNames_, i)
 263     {
 264         excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]);
 265
 266         if (excludedPatchLabels_[i] == -1)
 267         {
 268             FatalErrorIn
 269             (
 270                 "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
 271             )   << "cannot find MRF patch " << excludedPatchNames_[i]
 272                 << exit(FatalError);
 273         }
 274     }
 275
 276
 277     bool cellZoneFound = (cellZoneID_ != -1);
 278     reduce(cellZoneFound, orOp<bool>());
 279
 280     if (!cellZoneFound)
 281     {
 282         FatalErrorIn
 283         (
 284             "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
 285         )   << "cannot find MRF cellZone " << name_
 286             << exit(FatalError);
 287     }
 288
 289     setMRFFaces();
 290 }
 291
 292
 293 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 294
 295 void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
 296 {
 297     if (cellZoneID_ == -1)
 298     {
 299         return;
 300     }
 301
 302     const labelList& cells = mesh_.cellZones()[cellZoneID_];
 303     const scalarField& V = mesh_.V();
 304     vectorField& Usource = UEqn.source();
 305     const vectorField& U = UEqn.psi();
 306     const vector& Omega = Omega_.value();
 307
 308     forAll(cells, i)
 309     {
 310         label celli = cells[i];
 311         Usource[celli] -= V[celli]*(Omega ^ U[celli]);
 312     }
 313 }
 314
 315
 316 void Foam::MRFZone::addCoriolis
 317 (
 318     const volScalarField& rho,
 319     fvVectorMatrix& UEqn
 320 ) const
 321 {
 322     if (cellZoneID_ == -1)
 323     {
 324         return;
 325     }
 326
 327     const labelList& cells = mesh_.cellZones()[cellZoneID_];
 328     const scalarField& V = mesh_.V();
 329     vectorField& Usource = UEqn.source();
 330     const vectorField& U = UEqn.psi();
 331     const vector& Omega = Omega_.value();
 332
 333     forAll(cells, i)
 334     {
 335         label celli = cells[i];
 336         Usource[celli] -= V[celli]*rho[celli]*(Omega ^ U[celli]);
 337     }
 338 }
 339
 340
 341 void Foam::MRFZone::relativeVelocity(volVectorField& U) const
 342 {
 343     const volVectorField& C = mesh_.C();
 344
 345     const vector& origin = origin_.value();
 346     const vector& Omega = Omega_.value();
 347
 348     const labelList& cells = mesh_.cellZones()[cellZoneID_];
 349
 350     forAll(cells, i)
 351     {
 352         label celli = cells[i];
 353         U[celli] -= (Omega ^ (C[celli] - origin));
 354     }
 355
 356     // Included patches
 357     forAll(includedFaces_, patchi)
 358     {
 359         forAll(includedFaces_[patchi], i)
 360         {
 361             label patchFacei = includedFaces_[patchi][i];
 362             U.boundaryField()[patchi][patchFacei] = vector::zero;
 363         }
 364     }
 365
 366     // Excluded patches
 367     forAll(excludedFaces_, patchi)
 368     {
 369         forAll(excludedFaces_[patchi], i)
 370         {
 371             label patchFacei = excludedFaces_[patchi][i];
 372             U.boundaryField()[patchi][patchFacei] -=
 373                 (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
 374         }
 375     }
 376 }
 377
 378
 379 void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
 380 {
 381     relativeRhoFlux(geometricOneField(), phi);
 382 }
 383
 384
 385 void Foam::MRFZone::relativeFlux
 386 (
 387     const surfaceScalarField& rho,
 388     surfaceScalarField& phi
 389 ) const
 390 {
 391     relativeRhoFlux(rho, phi);
 392 }
 393
 394
 395 void Foam::MRFZone::absoluteFlux(surfaceScalarField& phi) const
 396 {
 397     absoluteRhoFlux(geometricOneField(), phi);
 398 }
 399
 400
 401 void Foam::MRFZone::absoluteFlux
 402 (
 403     const surfaceScalarField& rho,
 404     surfaceScalarField& phi
 405 ) const
 406 {
 407     absoluteRhoFlux(rho, phi);
 408 }
 409
 410
 411 void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
 412 {
 413     const vector& origin = origin_.value();
 414     const vector& Omega = Omega_.value();
 415
 416     // Included patches
 417     forAll(includedFaces_, patchi)
 418     {
 419         const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
 420
 421         vectorField pfld(U.boundaryField()[patchi]);
 422
 423         forAll(includedFaces_[patchi], i)
 424         {
 425             label patchFacei = includedFaces_[patchi][i];
 426
 427             pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin));
 428         }
 429
 430         U.boundaryField()[patchi] == pfld;
 431     }
 432 }
 433
 434
 435 // ************************************************************************* //