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) 2004-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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "MRFZone.H"
  27 #include "fvMesh.H"
  28 #include "volFields.H"
  29 #include "surfaceFields.H"
  30 #include "fvMatrices.H"
  31 #include "syncTools.H"
  32 #include "faceSet.H"
  33 #include "geometricOneField.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(Foam::MRFZone, 0);
  38
  39
  40 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  41
  42 void Foam::MRFZone::setMRFFaces()
  43 {
  44     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
  45
  46     // Type per face:
  47     //  0:not in zone
  48     //  1:moving with frame
  49     //  2:other
  50     labelList faceType(mesh_.nFaces(), 0);
  51
  52     // Determine faces in cell zone
  53     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  54     // (without constructing cells)
  55
  56     const labelList& own = mesh_.faceOwner();
  57     const labelList& nei = mesh_.faceNeighbour();
  58
  59     // Cells in zone
  60     boolList zoneCell(mesh_.nCells(), false);
  61
  62     if (cellZoneID_ != -1)
  63     {
  64         const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
  65         forAll(cellLabels, i)
  66         {
  67             zoneCell[cellLabels[i]] = true;
  68         }
  69     }
  70
  71
  72     label nZoneFaces = 0;
  73
  74     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
  75     {
  76         if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
  77         {
  78             faceType[faceI] = 1;
  79             nZoneFaces++;
  80         }
  81     }
  82
  83
  84     labelHashSet excludedPatches(excludedPatchLabels_);
  85
  86     forAll(patches, patchI)
  87     {
  88         const polyPatch& pp = patches[patchI];
  89
  90         if (pp.coupled() || excludedPatches.found(patchI))
  91         {
  92             forAll(pp, i)
  93             {
  94                 label faceI = pp.start()+i;
  95
  96                 if (zoneCell[own[faceI]])
  97                 {
  98                     faceType[faceI] = 2;
  99                     nZoneFaces++;
 100                 }
 101             }
 102         }
 103         else if (!isA<emptyPolyPatch>(pp))
 104         {
 105             forAll(pp, i)
 106             {
 107                 label faceI = pp.start()+i;
 108
 109                 if (zoneCell[own[faceI]])
 110                 {
 111                     faceType[faceI] = 1;
 112                     nZoneFaces++;
 113                 }
 114             }
 115         }
 116     }
 117
 118     // Now we have for faceType:
 119     //  0   : face not in cellZone
 120     //  1   : internal face or normal patch face
 121     //  2   : coupled patch face or excluded patch face
 122
 123     // Sort into lists per patch.
 124
 125     internalFaces_.setSize(mesh_.nFaces());
 126     label nInternal = 0;
 127
 128     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
 129     {
 130         if (faceType[faceI] == 1)
 131         {
 132             internalFaces_[nInternal++] = faceI;
 133         }
 134     }
 135     internalFaces_.setSize(nInternal);
 136
 137     labelList nIncludedFaces(patches.size(), 0);
 138     labelList nExcludedFaces(patches.size(), 0);
 139
 140     forAll(patches, patchi)
 141     {
 142         const polyPatch& pp = patches[patchi];
 143
 144         forAll(pp, patchFacei)
 145         {
 146             label faceI = pp.start()+patchFacei;
 147
 148             if (faceType[faceI] == 1)
 149             {
 150                 nIncludedFaces[patchi]++;
 151             }
 152             else if (faceType[faceI] == 2)
 153             {
 154                 nExcludedFaces[patchi]++;
 155             }
 156         }
 157     }
 158
 159     includedFaces_.setSize(patches.size());
 160     excludedFaces_.setSize(patches.size());
 161     forAll(nIncludedFaces, patchi)
 162     {
 163         includedFaces_[patchi].setSize(nIncludedFaces[patchi]);
 164         excludedFaces_[patchi].setSize(nExcludedFaces[patchi]);
 165     }
 166     nIncludedFaces = 0;
 167     nExcludedFaces = 0;
 168
 169     forAll(patches, patchi)
 170     {
 171         const polyPatch& pp = patches[patchi];
 172
 173         forAll(pp, patchFacei)
 174         {
 175             label faceI = pp.start()+patchFacei;
 176
 177             if (faceType[faceI] == 1)
 178             {
 179                 includedFaces_[patchi][nIncludedFaces[patchi]++] = patchFacei;
 180             }
 181             else if (faceType[faceI] == 2)
 182             {
 183                 excludedFaces_[patchi][nExcludedFaces[patchi]++] = patchFacei;
 184             }
 185         }
 186     }
 187
 188
 189     if (debug)
 190     {
 191         faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
 192         Pout<< "Writing " << internalFaces.size()
 193             << " internal faces in MRF zone to faceSet "
 194             << internalFaces.name() << endl;
 195         internalFaces.write();
 196
 197         faceSet MRFFaces(mesh_, "includedFaces", 100);
 198         forAll(includedFaces_, patchi)
 199         {
 200             forAll(includedFaces_[patchi], i)
 201             {
 202                 label patchFacei = includedFaces_[patchi][i];
 203                 MRFFaces.insert(patches[patchi].start()+patchFacei);
 204             }
 205         }
 206         Pout<< "Writing " << MRFFaces.size()
 207             << " patch faces in MRF zone to faceSet "
 208             << MRFFaces.name() << endl;
 209         MRFFaces.write();
 210
 211         faceSet excludedFaces(mesh_, "excludedFaces", 100);
 212         forAll(excludedFaces_, patchi)
 213         {
 214             forAll(excludedFaces_[patchi], i)
 215             {
 216                 label patchFacei = excludedFaces_[patchi][i];
 217                 excludedFaces.insert(patches[patchi].start()+patchFacei);
 218             }
 219         }
 220         Pout<< "Writing " << excludedFaces.size()
 221             << " faces in MRF zone with special handling to faceSet "
 222             << excludedFaces.name() << endl;
 223         excludedFaces.write();
 224     }
 225 }
 226
 227
 228 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 229
 230 Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
 231 :
 232     mesh_(mesh),
 233     name_(is),
 234     dict_(is),
 235     cellZoneID_(mesh_.cellZones().findZoneID(name_)),
 236     excludedPatchNames_
 237     (
 238         dict_.lookupOrDefault("nonRotatingPatches", wordList(0))
 239     ),
 240     origin_(dict_.lookup("origin")),
 241     axis_(dict_.lookup("axis")),
 242     omega_(dict_.lookup("omega")),
 243     Omega_("Omega", omega_*axis_)
 244 {
 245     if (dict_.found("patches"))
 246     {
 247         WarningIn("MRFZone(const fvMesh&, Istream&)")
 248             << "Ignoring entry 'patches'\n"
 249             << "    By default all patches within the rotating region rotate.\n"
 250             << "    Optionally supply excluded patches "
 251             << "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::absoluteVelocity(volVectorField& U) const
 380 {
 381     const volVectorField& C = mesh_.C();
 382
 383     const vector& origin = origin_.value();
 384     const vector& Omega = Omega_.value();
 385
 386     const labelList& cells = mesh_.cellZones()[cellZoneID_];
 387
 388     forAll(cells, i)
 389     {
 390         label celli = cells[i];
 391         U[celli] += (Omega ^ (C[celli] - origin));
 392     }
 393
 394     // Included patches
 395     forAll(includedFaces_, patchi)
 396     {
 397         forAll(includedFaces_[patchi], i)
 398         {
 399             label patchFacei = includedFaces_[patchi][i];
 400             U.boundaryField()[patchi][patchFacei] =
 401                 (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
 402         }
 403     }
 404
 405     // Excluded patches
 406     forAll(excludedFaces_, patchi)
 407     {
 408         forAll(excludedFaces_[patchi], i)
 409         {
 410             label patchFacei = excludedFaces_[patchi][i];
 411             U.boundaryField()[patchi][patchFacei] +=
 412                 (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
 413         }
 414     }
 415 }
 416
 417
 418 void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
 419 {
 420     relativeRhoFlux(geometricOneField(), phi);
 421 }
 422
 423
 424 void Foam::MRFZone::relativeFlux
 425 (
 426     const surfaceScalarField& rho,
 427     surfaceScalarField& phi
 428 ) const
 429 {
 430     relativeRhoFlux(rho, phi);
 431 }
 432
 433
 434 void Foam::MRFZone::absoluteFlux(surfaceScalarField& phi) const
 435 {
 436     absoluteRhoFlux(geometricOneField(), phi);
 437 }
 438
 439
 440 void Foam::MRFZone::absoluteFlux
 441 (
 442     const surfaceScalarField& rho,
 443     surfaceScalarField& phi
 444 ) const
 445 {
 446     absoluteRhoFlux(rho, phi);
 447 }
 448
 449
 450 void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
 451 {
 452     const vector& origin = origin_.value();
 453     const vector& Omega = Omega_.value();
 454
 455     // Included patches
 456     forAll(includedFaces_, patchi)
 457     {
 458         const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
 459
 460         vectorField pfld(U.boundaryField()[patchi]);
 461
 462         forAll(includedFaces_[patchi], i)
 463         {
 464             label patchFacei = includedFaces_[patchi][i];
 465
 466             pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin));
 467         }
 468
 469         U.boundaryField()[patchi] == pfld;
 470     }
 471 }
 472
 473
 474 // ************************************************************************* //