src/finiteVolume/cfdTools/general/porousMedia/porousZone.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-2011 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 "porousZone.H"
  27 #include "fvMesh.H"
  28 #include "fvMatrices.H"
  29 #include "geometricOneField.H"
  30
  31 // * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
  32
  33 // adjust negative resistance values to be multiplier of max value
  34 void Foam::porousZone::adjustNegativeResistance(dimensionedVector& resist)
  35 {
  36     scalar maxCmpt = max(0, cmptMax(resist.value()));
  37
  38     if (maxCmpt < 0)
  39     {
  40         FatalErrorIn
  41         (
  42             "Foam::porousZone::porousZone::adjustNegativeResistance"
  43             "(dimensionedVector&)"
  44         )   << "negative resistances! " << resist
  45             << exit(FatalError);
  46     }
  47     else
  48     {
  49         vector& val = resist.value();
  50         for (label cmpt=0; cmpt < vector::nComponents; ++cmpt)
  51         {
  52             if (val[cmpt] < 0)
  53             {
  54                 val[cmpt] *= -maxCmpt;
  55             }
  56         }
  57     }
  58 }
  59
  60
  61 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  62
  63 Foam::porousZone::porousZone
  64 (
  65     const keyType& key,
  66     const fvMesh& mesh,
  67     const dictionary& dict
  68 )
  69 :
  70     key_(key),
  71     mesh_(mesh),
  72     dict_(dict),
  73     cellZoneIds_(mesh_.cellZones().findIndices(key_)),
  74     coordSys_(dict, mesh),
  75     porosity_(1),
  76     intensity_(0),
  77     mixingLength_(0),
  78     C0_(0),
  79     C1_(0),
  80     D_("D", dimensionSet(0, -2, 0, 0, 0), tensor::zero),
  81     F_("F", dimensionSet(0, -1, 0, 0, 0), tensor::zero)
  82 {
  83     Info<< "Creating porous zone: " << key_ << endl;
  84
  85     bool foundZone = !cellZoneIds_.empty();
  86     reduce(foundZone, orOp<bool>());
  87
  88     if (!foundZone && Pstream::master())
  89     {
  90         FatalErrorIn
  91         (
  92             "Foam::porousZone::porousZone"
  93             "(const keyType&, const fvMesh&, const dictionary&)"
  94         )   << "cannot find porous cellZone " << key_
  95             << exit(FatalError);
  96     }
  97
  98
  99     // porosity
 100     if
 101     (
 102         dict_.readIfPresent("porosity", porosity_)
 103      && (porosity_ <= 0.0 || porosity_ > 1.0)
 104     )
 105     {
 106         FatalIOErrorIn
 107         (
 108             "Foam::porousZone::porousZone"
 109             "(const keyType&, const fvMesh&, const dictionary&)",
 110             dict_
 111         )
 112             << "out-of-range porosity value " << porosity_
 113             << exit(FatalIOError);
 114     }
 115
 116     // turbulent intensity
 117     if
 118     (
 119         dict_.readIfPresent("intensity", intensity_)
 120      && (intensity_ <= 0.0 || intensity_ > 1.0)
 121     )
 122     {
 123         FatalIOErrorIn
 124         (
 125             "Foam::porousZone::porousZone"
 126             "(const keyType&, const fvMesh&, const dictionary&)",
 127             dict_
 128         )
 129             << "out-of-range turbulent intensity value " << intensity_
 130             << exit(FatalIOError);
 131     }
 132
 133     // turbulent length scale
 134     if
 135     (
 136         dict_.readIfPresent("mixingLength", mixingLength_)
 137      && (mixingLength_ <= 0.0)
 138     )
 139     {
 140         FatalIOErrorIn
 141         (
 142             "Foam::porousZone::porousZone"
 143             "(const keyType&, const fvMesh&, const dictionary&)",
 144             dict_
 145         )
 146             << "out-of-range turbulent length scale " << mixingLength_
 147             << exit(FatalIOError);
 148     }
 149
 150
 151     // powerLaw coefficients
 152     if (const dictionary* dictPtr = dict_.subDictPtr("powerLaw"))
 153     {
 154         dictPtr->readIfPresent("C0", C0_);
 155         dictPtr->readIfPresent("C1", C1_);
 156     }
 157
 158     // Darcy-Forchheimer coefficients
 159     if (const dictionary* dictPtr = dict_.subDictPtr("Darcy"))
 160     {
 161         // local-to-global transformation tensor
 162         const tensor& E = coordSys_.R();
 163
 164         dimensionedVector d(vector::zero);
 165         if (dictPtr->readIfPresent("d", d))
 166         {
 167             if (D_.dimensions() != d.dimensions())
 168             {
 169                 FatalIOErrorIn
 170                 (
 171                     "Foam::porousZone::porousZone"
 172                     "(const keyType&, const fvMesh&, const dictionary&)",
 173                     dict_
 174                 )   << "incorrect dimensions for d: " << d.dimensions()
 175                     << " should be " << D_.dimensions()
 176                     << exit(FatalIOError);
 177             }
 178
 179             adjustNegativeResistance(d);
 180
 181             D_.value().xx() = d.value().x();
 182             D_.value().yy() = d.value().y();
 183             D_.value().zz() = d.value().z();
 184             D_.value() = (E & D_ & E.T()).value();
 185         }
 186
 187         dimensionedVector f(vector::zero);
 188         if (dictPtr->readIfPresent("f", f))
 189         {
 190             if (F_.dimensions() != f.dimensions())
 191             {
 192                 FatalIOErrorIn
 193                 (
 194                     "Foam::porousZone::porousZone"
 195                     "(const keyType&, const fvMesh&, const dictionary&)",
 196                     dict_
 197                 )   << "incorrect dimensions for f: " << f.dimensions()
 198                     << " should be " << F_.dimensions()
 199                     << exit(FatalIOError);
 200             }
 201
 202             adjustNegativeResistance(f);
 203
 204             // leading 0.5 is from 1/2 * rho
 205             F_.value().xx() = 0.5*f.value().x();
 206             F_.value().yy() = 0.5*f.value().y();
 207             F_.value().zz() = 0.5*f.value().z();
 208             F_.value() = (E & F_ & E.T()).value();
 209         }
 210     }
 211
 212     // it is an error not to define anything
 213     if
 214     (
 215         C0_ <= VSMALL
 216      && magSqr(D_.value()) <= VSMALL
 217      && magSqr(F_.value()) <= VSMALL
 218     )
 219     {
 220         FatalIOErrorIn
 221         (
 222             "Foam::porousZone::porousZone"
 223             "(const keyType&, const fvMesh&, const dictionary&)",
 224             dict_
 225         )   << "neither powerLaw (C0/C1) "
 226                "nor Darcy-Forchheimer law (d/f) specified"
 227             << exit(FatalIOError);
 228     }
 229
 230     // feedback for the user
 231     if (dict.lookupOrDefault("printCoeffs", false))
 232     {
 233         writeDict(Info, false);
 234     }
 235 }
 236
 237
 238 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 239
 240 void Foam::porousZone::addResistance(fvVectorMatrix& UEqn) const
 241 {
 242     if (cellZoneIds_.empty())
 243     {
 244         return;
 245     }
 246
 247     bool compressible = false;
 248     if (UEqn.dimensions() == dimensionSet(1, 1, -2, 0, 0))
 249     {
 250         compressible = true;
 251     }
 252
 253     const scalarField& V = mesh_.V();
 254     scalarField& Udiag = UEqn.diag();
 255     vectorField& Usource = UEqn.source();
 256     const vectorField& U = UEqn.psi();
 257
 258     if (C0_ > VSMALL)
 259     {
 260         if (compressible)
 261         {
 262             addPowerLawResistance
 263             (
 264                 Udiag,
 265                 V,
 266                 mesh_.lookupObject<volScalarField>("rho"),
 267                 U
 268             );
 269         }
 270         else
 271         {
 272             addPowerLawResistance
 273             (
 274                 Udiag,
 275                 V,
 276                 geometricOneField(),
 277                 U
 278             );
 279         }
 280     }
 281
 282     const tensor& D = D_.value();
 283     const tensor& F = F_.value();
 284
 285     if (magSqr(D) > VSMALL || magSqr(F) > VSMALL)
 286     {
 287         if (compressible)
 288         {
 289             addViscousInertialResistance
 290             (
 291                 Udiag,
 292                 Usource,
 293                 V,
 294                 mesh_.lookupObject<volScalarField>("rho"),
 295                 mesh_.lookupObject<volScalarField>("mu"),
 296                 U
 297             );
 298         }
 299         else
 300         {
 301             addViscousInertialResistance
 302             (
 303                 Udiag,
 304                 Usource,
 305                 V,
 306                 geometricOneField(),
 307                 mesh_.lookupObject<volScalarField>("nu"),
 308                 U
 309             );
 310         }
 311     }
 312 }
 313
 314
 315 void Foam::porousZone::addResistance
 316 (
 317     const fvVectorMatrix& UEqn,
 318     volTensorField& AU,
 319     bool correctAUprocBC
 320 ) const
 321 {
 322     if (cellZoneIds_.empty())
 323     {
 324         return;
 325     }
 326
 327     bool compressible = false;
 328     if (UEqn.dimensions() == dimensionSet(1, 1, -2, 0, 0))
 329     {
 330         compressible = true;
 331     }
 332
 333     const vectorField& U = UEqn.psi();
 334
 335     if (C0_ > VSMALL)
 336     {
 337         if (compressible)
 338         {
 339             addPowerLawResistance
 340             (
 341                 AU,
 342                 mesh_.lookupObject<volScalarField>("rho"),
 343                 U
 344             );
 345         }
 346         else
 347         {
 348             addPowerLawResistance
 349             (
 350                 AU,
 351                 geometricOneField(),
 352                 U
 353             );
 354         }
 355     }
 356
 357     const tensor& D = D_.value();
 358     const tensor& F = F_.value();
 359
 360     if (magSqr(D) > VSMALL || magSqr(F) > VSMALL)
 361     {
 362         if (compressible)
 363         {
 364             addViscousInertialResistance
 365             (
 366                 AU,
 367                 mesh_.lookupObject<volScalarField>("rho"),
 368                 mesh_.lookupObject<volScalarField>("mu"),
 369                 U
 370             );
 371         }
 372         else
 373         {
 374             addViscousInertialResistance
 375             (
 376                 AU,
 377                 geometricOneField(),
 378                 mesh_.lookupObject<volScalarField>("nu"),
 379                 U
 380             );
 381         }
 382     }
 383
 384     if (correctAUprocBC)
 385     {
 386         // Correct the boundary conditions of the tensorial diagonal to ensure
 387         // processor boundaries are correctly handled when AU^-1 is interpolated
 388         // for the pressure equation.
 389         AU.correctBoundaryConditions();
 390     }
 391 }
 392
 393
 394 void Foam::porousZone::writeDict(Ostream& os, bool subDict) const
 395 {
 396     if (subDict)
 397     {
 398         os  << indent << token::BEGIN_BLOCK << incrIndent << nl;
 399         os.writeKeyword("name")
 400             << zoneName() << token::END_STATEMENT << nl;
 401     }
 402     else
 403     {
 404         os  << indent << zoneName() << nl
 405             << indent << token::BEGIN_BLOCK << incrIndent << nl;
 406     }
 407
 408     if (dict_.found("note"))
 409     {
 410         os.writeKeyword("note")
 411             << string(dict_.lookup("note")) << token::END_STATEMENT << nl;
 412     }
 413
 414     coordSys_.writeDict(os, true);
 415
 416     if (dict_.found("porosity"))
 417     {
 418         os.writeKeyword("porosity")
 419             << porosity() << token::END_STATEMENT << nl;
 420     }
 421
 422     if (dict_.found("intensity"))
 423     {
 424         os.writeKeyword("intensity")
 425             << intensity() << token::END_STATEMENT << nl;
 426     }
 427
 428     if (dict_.found("mixingLength"))
 429     {
 430         os.writeKeyword("mixingLength")
 431             << mixingLength() << token::END_STATEMENT << nl;
 432     }
 433
 434     // powerLaw coefficients
 435     if (const dictionary* dictPtr = dict_.subDictPtr("powerLaw"))
 436     {
 437         os  << indent << "powerLaw";
 438         dictPtr->write(os);
 439     }
 440
 441     // Darcy-Forchheimer coefficients
 442     if (const dictionary* dictPtr = dict_.subDictPtr("Darcy"))
 443     {
 444         os  << indent << "Darcy";
 445         dictPtr->write(os);
 446     }
 447
 448     // thermalModel
 449     if (const dictionary* dictPtr = dict_.subDictPtr("thermalModel"))
 450     {
 451         os  << indent << "thermalModel";
 452         dictPtr->write(os);
 453     }
 454
 455     os  << decrIndent << indent << token::END_BLOCK << endl;
 456 }
 457
 458
 459 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 460
 461 Foam::Ostream& Foam::operator<<(Ostream& os, const porousZone& pz)
 462 {
 463     pz.writeDict(os);
 464     return os;
 465 }
 466
 467 // ************************************************************************* //