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