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

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