applications/utilities/postProcessing/velocityField/Pe/Pe.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Application
  25     Pe
  26
  27 Description
  28     Calculates and writes the Pe number as a surfaceScalarField obtained from
  29     field phi.
  30
  31     The -noWrite option just outputs the max/min values without writing
  32     the field.
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "calc.H"
  37 #include "fvc.H"
  38
  39 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
  40 #include "incompressible/RAS/RASModel/RASModel.H"
  41 #include "incompressible/LES/LESModel/LESModel.H"
  42 #include "basicPsiThermo.H"
  43 #include "compressible/RAS/RASModel/RASModel.H"
  44 #include "compressible/LES/LESModel/LESModel.H"
  45
  46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  47
  48 void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
  49 {
  50     bool writeResults = !args.optionFound("noWrite");
  51
  52     IOobject phiHeader
  53     (
  54         "phi",
  55         runTime.timeName(),
  56         mesh,
  57         IOobject::MUST_READ
  58     );
  59
  60     if (phiHeader.headerOk())
  61     {
  62         autoPtr<surfaceScalarField> PePtr;
  63
  64         Info<< "    Reading phi" << endl;
  65         surfaceScalarField phi(phiHeader, mesh);
  66
  67         volVectorField U
  68         (
  69             IOobject
  70             (
  71                 "U",
  72                 runTime.timeName(),
  73                 mesh,
  74                 IOobject::MUST_READ
  75             ),
  76             mesh
  77         );
  78
  79         IOobject RASPropertiesHeader
  80         (
  81             "RASProperties",
  82             runTime.constant(),
  83             mesh,
  84             IOobject::MUST_READ,
  85             IOobject::NO_WRITE
  86         );
  87
  88         IOobject LESPropertiesHeader
  89         (
  90             "LESProperties",
  91             runTime.constant(),
  92             mesh,
  93             IOobject::MUST_READ,
  94             IOobject::NO_WRITE
  95         );
  96
  97         Info<< "    Calculating Pe" << endl;
  98
  99         if (phi.dimensions() == dimensionSet(0, 3, -1, 0, 0))
 100         {
 101             if (RASPropertiesHeader.headerOk())
 102             {
 103                 IOdictionary RASProperties(RASPropertiesHeader);
 104
 105                 singlePhaseTransportModel laminarTransport(U, phi);
 106
 107                 autoPtr<incompressible::RASModel> RASModel
 108                 (
 109                     incompressible::RASModel::New
 110                     (
 111                         U,
 112                         phi,
 113                         laminarTransport
 114                     )
 115                 );
 116
 117                 PePtr.set
 118                 (
 119                     new surfaceScalarField
 120                     (
 121                         IOobject
 122                         (
 123                             "Pe",
 124                             runTime.timeName(),
 125                             mesh,
 126                             IOobject::NO_READ
 127                         ),
 128                         mag(phi)
 129                        /(
 130                             mesh.magSf()
 131                           * mesh.surfaceInterpolation::deltaCoeffs()
 132                           * fvc::interpolate(RASModel->nuEff())
 133                          )
 134                     )
 135                 );
 136             }
 137             else if (LESPropertiesHeader.headerOk())
 138             {
 139                 IOdictionary LESProperties(LESPropertiesHeader);
 140
 141                 singlePhaseTransportModel laminarTransport(U, phi);
 142
 143                 autoPtr<incompressible::LESModel> sgsModel
 144                 (
 145                     incompressible::LESModel::New(U, phi, laminarTransport)
 146                 );
 147
 148                 PePtr.set
 149                 (
 150                     new surfaceScalarField
 151                     (
 152                         IOobject
 153                         (
 154                             "Pe",
 155                             runTime.timeName(),
 156                             mesh,
 157                             IOobject::NO_READ
 158                         ),
 159                         mag(phi)
 160                        /(
 161                             mesh.magSf()
 162                           * mesh.surfaceInterpolation::deltaCoeffs()
 163                           * fvc::interpolate(sgsModel->nuEff())
 164                         )
 165                     )
 166                 );
 167             }
 168             else
 169             {
 170                 IOdictionary transportProperties
 171                 (
 172                     IOobject
 173                     (
 174                         "transportProperties",
 175                         runTime.constant(),
 176                         mesh,
 177                         IOobject::MUST_READ,
 178                         IOobject::NO_WRITE
 179                     )
 180                 );
 181
 182                 dimensionedScalar nu(transportProperties.lookup("nu"));
 183
 184                 PePtr.set
 185                 (
 186                     new surfaceScalarField
 187                     (
 188                         IOobject
 189                         (
 190                             "Pe",
 191                             runTime.timeName(),
 192                             mesh,
 193                             IOobject::NO_READ
 194                         ),
 195                         mesh.surfaceInterpolation::deltaCoeffs()
 196                       * (mag(phi)/mesh.magSf())*(runTime.deltaT()/nu)
 197                     )
 198                 );
 199             }
 200         }
 201         else if (phi.dimensions() == dimensionSet(1, 0, -1, 0, 0))
 202         {
 203             if (RASPropertiesHeader.headerOk())
 204             {
 205                 IOdictionary RASProperties(RASPropertiesHeader);
 206
 207                 autoPtr<basicPsiThermo> thermo(basicPsiThermo::New(mesh));
 208
 209                 volScalarField rho
 210                 (
 211                     IOobject
 212                     (
 213                         "rho",
 214                         runTime.timeName(),
 215                         mesh
 216                     ),
 217                     thermo->rho()
 218                 );
 219
 220                 autoPtr<compressible::RASModel> RASModel
 221                 (
 222                     compressible::RASModel::New
 223                     (
 224                         rho,
 225                         U,
 226                         phi,
 227                         thermo()
 228                     )
 229                 );
 230
 231                 PePtr.set
 232                 (
 233                     new surfaceScalarField
 234                     (
 235                         IOobject
 236                         (
 237                             "Pe",
 238                             runTime.timeName(),
 239                             mesh,
 240                             IOobject::NO_READ
 241                         ),
 242                         mag(phi)
 243                        /(
 244                             mesh.magSf()
 245                           * mesh.surfaceInterpolation::deltaCoeffs()
 246                           * fvc::interpolate(RASModel->muEff())
 247                         )
 248                     )
 249                 );
 250             }
 251             else if (LESPropertiesHeader.headerOk())
 252             {
 253                 IOdictionary LESProperties(LESPropertiesHeader);
 254
 255                 autoPtr<basicPsiThermo> thermo(basicPsiThermo::New(mesh));
 256
 257                 volScalarField rho
 258                 (
 259                     IOobject
 260                     (
 261                         "rho",
 262                         runTime.timeName(),
 263                         mesh
 264                     ),
 265                     thermo->rho()
 266                 );
 267
 268                 autoPtr<compressible::LESModel> sgsModel
 269                 (
 270                     compressible::LESModel::New(rho, U, phi, thermo())
 271                 );
 272
 273                 PePtr.set
 274                 (
 275                     new surfaceScalarField
 276                     (
 277                         IOobject
 278                         (
 279                             "Pe",
 280                             runTime.timeName(),
 281                             mesh,
 282                             IOobject::NO_READ
 283                         ),
 284                         mag(phi)
 285                        /(
 286                             mesh.magSf()
 287                           * mesh.surfaceInterpolation::deltaCoeffs()
 288                           * fvc::interpolate(sgsModel->muEff())
 289                         )
 290                     )
 291                 );
 292             }
 293             else
 294             {
 295                 IOdictionary transportProperties
 296                 (
 297                     IOobject
 298                     (
 299                         "transportProperties",
 300                         runTime.constant(),
 301                         mesh,
 302                         IOobject::MUST_READ,
 303                         IOobject::NO_WRITE
 304                     )
 305                 );
 306
 307                 dimensionedScalar mu(transportProperties.lookup("mu"));
 308
 309                 PePtr.set
 310                 (
 311                     new surfaceScalarField
 312                     (
 313                         IOobject
 314                         (
 315                             "Pe",
 316                             runTime.timeName(),
 317                             mesh,
 318                             IOobject::NO_READ
 319                         ),
 320                         mesh.surfaceInterpolation::deltaCoeffs()
 321                       * (mag(phi)/(mesh.magSf()))*(runTime.deltaT()/mu)
 322                     )
 323                 );
 324             }
 325         }
 326         else
 327         {
 328             FatalErrorIn(args.executable())
 329                 << "Incorrect dimensions of phi: " << phi.dimensions()
 330                     << abort(FatalError);
 331         }
 332
 333
 334         // can also check how many cells exceed a particular Pe limit
 335         /*
 336         {
 337             label count = 0;
 338             label PeLimit = 200;
 339             forAll(PePtr(), i)
 340             {
 341                 if (PePtr()[i] > PeLimit)
 342                 {
 343                     count++;
 344                 }
 345
 346             }
 347
 348             Info<< "Fraction > " << PeLimit << " = "
 349                 << scalar(count)/Pe.size() << endl;
 350         }
 351         */
 352
 353         Info << "Pe max : " << max(PePtr()).value() << endl;
 354
 355         if (writeResults)
 356         {
 357             PePtr().write();
 358         }
 359     }
 360     else
 361     {
 362         Info<< "    No phi" << endl;
 363     }
 364
 365     Info<< "\nEnd\n" << endl;
 366 }
 367
 368
 369 // ************************************************************************* //