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-2009 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
  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 using 'nonRotatingPatches'."
 251             << endl;
 252     }
 253
 254     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 255
 256     axis_ = axis_/mag(axis_);
 257     Omega_ = omega_*axis_;
 258
 259     excludedPatchLabels_.setSize(excludedPatchNames_.size());
 260
 261     forAll(excludedPatchNames_, i)
 262     {
 263         excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]);
 264
 265         if (excludedPatchLabels_[i] == -1)
 266         {
 267             FatalErrorIn
 268             (
 269                 "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
 270             )   << "cannot find MRF patch " << excludedPatchNames_[i]
 271                 << exit(FatalError);
 272         }
 273     }
 274
 275
 276     bool cellZoneFound = (cellZoneID_ != -1);
 277     reduce(cellZoneFound, orOp<bool>());
 278
 279     if (!cellZoneFound)
 280     {
 281         FatalErrorIn
 282         (
 283             "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
 284         )   << "cannot find MRF cellZone " << name_
 285             << exit(FatalError);
 286     }
 287
 288     setMRFFaces();
 289 }
 290
 291
 292 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 293
 294 void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
 295 {
 296     if (cellZoneID_ == -1)
 297     {
 298         return;
 299     }
 300
 301     const labelList& cells = mesh_.cellZones()[cellZoneID_];
 302     const scalarField& V = mesh_.V();
 303     vectorField& Usource = UEqn.source();
 304     const vectorField& U = UEqn.psi();
 305     const vector& Omega = Omega_.value();
 306
 307     forAll(cells, i)
 308     {
 309         label celli = cells[i];
 310         Usource[celli] -= V[celli]*(Omega ^ U[celli]);
 311     }
 312 }
 313
 314
 315 void Foam::MRFZone::addCoriolis
 316 (
 317     const volScalarField& rho,
 318     fvVectorMatrix& UEqn
 319 ) const
 320 {
 321     if (cellZoneID_ == -1)
 322     {
 323         return;
 324     }
 325
 326     const labelList& cells = mesh_.cellZones()[cellZoneID_];
 327     const scalarField& V = mesh_.V();
 328     vectorField& Usource = UEqn.source();
 329     const vectorField& U = UEqn.psi();
 330     const vector& Omega = Omega_.value();
 331
 332     forAll(cells, i)
 333     {
 334         label celli = cells[i];
 335         Usource[celli] -= V[celli]*rho[celli]*(Omega ^ U[celli]);
 336     }
 337 }
 338
 339
 340 void Foam::MRFZone::relativeVelocity(volVectorField& U) const
 341 {
 342     const volVectorField& C = mesh_.C();
 343
 344     const vector& origin = origin_.value();
 345     const vector& Omega = Omega_.value();
 346
 347     const labelList& cells = mesh_.cellZones()[cellZoneID_];
 348
 349     forAll(cells, i)
 350     {
 351         label celli = cells[i];
 352         U[celli] -= (Omega ^ (C[celli] - origin));
 353     }
 354
 355     // Included patches
 356     forAll(includedFaces_, patchi)
 357     {
 358         forAll(includedFaces_[patchi], i)
 359         {
 360             label patchFacei = includedFaces_[patchi][i];
 361             U.boundaryField()[patchi][patchFacei] = vector::zero;
 362         }
 363     }
 364
 365     // Excluded patches
 366     forAll(excludedFaces_, patchi)
 367     {
 368         forAll(excludedFaces_[patchi], i)
 369         {
 370             label patchFacei = excludedFaces_[patchi][i];
 371             U.boundaryField()[patchi][patchFacei] -=
 372                 (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
 373         }
 374     }
 375 }
 376
 377
 378 void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
 379 {
 380     const surfaceVectorField& Cf = mesh_.Cf();
 381     const surfaceVectorField& Sf = mesh_.Sf();
 382
 383     const vector& origin = origin_.value();
 384     const vector& Omega = Omega_.value();
 385
 386     // Internal faces
 387     forAll(internalFaces_, i)
 388     {
 389         label facei = internalFaces_[i];
 390         phi[facei] -= (Omega ^ (Cf[facei] - origin)) & Sf[facei];
 391     }
 392
 393     // Included patches
 394     forAll(includedFaces_, patchi)
 395     {
 396         forAll(includedFaces_[patchi], i)
 397         {
 398             label patchFacei = includedFaces_[patchi][i];
 399
 400             phi.boundaryField()[patchi][patchFacei] = 0.0;
 401         }
 402     }
 403
 404     // Excluded patches
 405     forAll(excludedFaces_, patchi)
 406     {
 407         forAll(excludedFaces_[patchi], i)
 408         {
 409             label patchFacei = excludedFaces_[patchi][i];
 410
 411             phi.boundaryField()[patchi][patchFacei] -=
 412                 (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
 413               & Sf.boundaryField()[patchi][patchFacei];
 414         }
 415     }
 416 }
 417
 418
 419 void Foam::MRFZone::absoluteFlux(surfaceScalarField& phi) const
 420 {
 421     const surfaceVectorField& Cf = mesh_.Cf();
 422     const surfaceVectorField& Sf = mesh_.Sf();
 423
 424     const vector& origin = origin_.value();
 425     const vector& Omega = Omega_.value();
 426
 427     // Internal faces
 428     forAll(internalFaces_, i)
 429     {
 430         label facei = internalFaces_[i];
 431         phi[facei] += (Omega ^ (Cf[facei] - origin)) & Sf[facei];
 432     }
 433
 434     // Included patches
 435     forAll(includedFaces_, patchi)
 436     {
 437         forAll(includedFaces_[patchi], i)
 438         {
 439             label patchFacei = includedFaces_[patchi][i];
 440
 441             phi.boundaryField()[patchi][patchFacei] +=
 442                 (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
 443               & Sf.boundaryField()[patchi][patchFacei];
 444         }
 445     }
 446
 447     // Excluded patches
 448     forAll(excludedFaces_, patchi)
 449     {
 450         forAll(excludedFaces_[patchi], i)
 451         {
 452             label patchFacei = excludedFaces_[patchi][i];
 453
 454             phi.boundaryField()[patchi][patchFacei] +=
 455                 (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
 456               & Sf.boundaryField()[patchi][patchFacei];
 457         }
 458     }
 459 }
 460
 461
 462 void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
 463 {
 464     const vector& origin = origin_.value();
 465     const vector& Omega = Omega_.value();
 466
 467     // Included patches
 468     forAll(includedFaces_, patchi)
 469     {
 470         const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
 471
 472         vectorField pfld(U.boundaryField()[patchi]);
 473
 474         forAll(includedFaces_[patchi], i)
 475         {
 476             label patchFacei = includedFaces_[patchi][i];
 477
 478             pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin));
 479         }
 480
 481         U.boundaryField()[patchi] == pfld;
 482     }
 483 }
 484
 485
 486 // ************************************************************************* //