applications/utilities/postProcessing/wall/yPlusRAS/yPlusRAS.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     yPlusRAS
  26
  27 Description
  28     Calculates and reports yPlus for all wall patches, for the specified times
  29     when using RAS turbulence models.
  30
  31     Default behaviour assumes operating in incompressible mode. To apply to
  32     compressible RAS cases, use the -compressible option.
  33
  34     Extended version for being able to handle two phase flows using the
  35     -twoPhase option.
  36
  37     Frank Albina, 16/Nov/2009
  38
  39 \*---------------------------------------------------------------------------*/
  40
  41 #include "fvCFD.H"
  42 #include "incompressible/incompressibleTwoPhaseMixture/twoPhaseMixture.H"
  43 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
  44 #include "incompressible/RAS/RASModel/RASModel.H"
  45 #include "nutWallFunction/nutWallFunctionFvPatchScalarField.H"
  46
  47 #include "basicPsiThermo.H"
  48 #include "compressible/RAS/RASModel/RASModel.H"
  49 #include "mutWallFunction/mutWallFunctionFvPatchScalarField.H"
  50
  51 #include "wallDist.H"
  52
  53 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  54
  55 void calcIncompressibleYPlus
  56 (
  57     const fvMesh& mesh,
  58     const Time& runTime,
  59     const volVectorField& U,
  60     volScalarField& yPlus
  61 )
  62 {
  63     typedef incompressible::RASModels::nutWallFunctionFvPatchScalarField
  64         wallFunctionPatchField;
  65
  66 #   include "createPhi.H"
  67
  68     singlePhaseTransportModel laminarTransport(U, phi);
  69
  70     autoPtr<incompressible::RASModel> RASModel
  71     (
  72         incompressible::RASModel::New(U, phi, laminarTransport)
  73     );
  74
  75     const volScalarField::GeometricBoundaryField nutPatches =
  76         RASModel->nut()().boundaryField();
  77
  78     bool foundNutPatch = false;
  79     forAll(nutPatches, patchi)
  80     {
  81         if (isA<wallFunctionPatchField>(nutPatches[patchi]))
  82         {
  83             foundNutPatch = true;
  84
  85             const wallFunctionPatchField& nutPw =
  86                 dynamic_cast<const wallFunctionPatchField&>
  87                     (nutPatches[patchi]);
  88
  89             yPlus.boundaryField()[patchi] = nutPw.yPlus();
  90             const scalarField& Yp = yPlus.boundaryField()[patchi];
  91
  92             Info<< "Patch " << patchi
  93                 << " named " << nutPw.patch().name()
  94                 << " y+ : min: " << min(Yp) << " max: " << max(Yp)
  95                 << " average: " << average(Yp) << nl << endl;
  96         }
  97     }
  98
  99     if (!foundNutPatch)
 100     {
 101         Info<< "    no " << wallFunctionPatchField::typeName << " patches"
 102             << endl;
 103     }
 104 }
 105
 106
 107 void calcCompressibleYPlus
 108 (
 109     const fvMesh& mesh,
 110     const Time& runTime,
 111     const volVectorField& U,
 112     volScalarField& yPlus
 113 )
 114 {
 115     typedef compressible::RASModels::mutWallFunctionFvPatchScalarField
 116         wallFunctionPatchField;
 117
 118     IOobject rhoHeader
 119     (
 120         "rho",
 121         runTime.timeName(),
 122         mesh,
 123         IOobject::MUST_READ,
 124         IOobject::NO_WRITE
 125     );
 126
 127     if (!rhoHeader.headerOk())
 128     {
 129         Info<< "    no rho field" << endl;
 130         return;
 131     }
 132
 133     Info << "Reading field rho\n" << endl;
 134     volScalarField rho(rhoHeader, mesh);
 135
 136 #   include "compressibleCreatePhi.H"
 137
 138     autoPtr<basicPsiThermo> pThermo
 139     (
 140         basicPsiThermo::New(mesh)
 141     );
 142     basicPsiThermo& thermo = pThermo();
 143
 144     autoPtr<compressible::RASModel> RASModel
 145     (
 146         compressible::RASModel::New
 147         (
 148             rho,
 149             U,
 150             phi,
 151             thermo
 152         )
 153     );
 154
 155     const volScalarField::GeometricBoundaryField mutPatches =
 156         RASModel->mut()().boundaryField();
 157
 158     bool foundMutPatch = false;
 159     forAll(mutPatches, patchi)
 160     {
 161         if (isA<wallFunctionPatchField>(mutPatches[patchi]))
 162         {
 163             foundMutPatch = true;
 164
 165             const wallFunctionPatchField& mutPw =
 166                 dynamic_cast<const wallFunctionPatchField&>
 167                     (mutPatches[patchi]);
 168
 169             yPlus.boundaryField()[patchi] = mutPw.yPlus();
 170             const scalarField& Yp = yPlus.boundaryField()[patchi];
 171
 172             Info<< "Patch " << patchi
 173                 << " named " << mutPw.patch().name()
 174                 << " y+ : min: " << min(Yp) << " max: " << max(Yp)
 175                 << " average: " << average(Yp) << nl << endl;
 176         }
 177     }
 178
 179     if (!foundMutPatch)
 180     {
 181         Info<< "    no " << wallFunctionPatchField::typeName << " patches"
 182             << endl;
 183     }
 184 }
 185
 186
 187 // Calculate two phase Y+
 188 void calcTwoPhaseYPlus
 189 (
 190     const fvMesh& mesh,
 191     const Time& runTime,
 192     const volVectorField& U,
 193     volScalarField& yPlus
 194 )
 195 {
 196     typedef incompressible::RASModels::nutWallFunctionFvPatchScalarField
 197         wallFunctionPatchField;
 198
 199 #   include "createPhi.H"
 200
 201     IOobject alphaHeader
 202     (
 203         "alpha1",
 204         runTime.timeName(),
 205         mesh,
 206         IOobject::MUST_READ,
 207         IOobject::NO_WRITE
 208     );
 209
 210     if (!alphaHeader.headerOk())
 211     {
 212         Info<< "    no alpha1 field" << endl;
 213         return;
 214     }
 215
 216     Info << "Reading field alpha1\n" << endl;
 217     volScalarField alpha(alphaHeader, mesh);
 218
 219     Info<< "Reading transportProperties\n" << endl;
 220     twoPhaseMixture twoPhaseProperties(U, phi, "alpha1");
 221
 222     autoPtr<incompressible::RASModel> RASModel
 223     (
 224         incompressible::RASModel::New(U, phi, twoPhaseProperties)
 225     );
 226
 227     const volScalarField::GeometricBoundaryField nutPatches =
 228         RASModel->nut()().boundaryField();
 229
 230     bool foundNutPatch = false;
 231     forAll(nutPatches, patchi)
 232     {
 233         if (isA<wallFunctionPatchField>(nutPatches[patchi]))
 234         {
 235             foundNutPatch = true;
 236
 237             const wallFunctionPatchField& nutPw =
 238                 dynamic_cast<const wallFunctionPatchField&>
 239                     (nutPatches[patchi]);
 240
 241             yPlus.boundaryField()[patchi] = nutPw.yPlus();
 242             const scalarField& Yp = yPlus.boundaryField()[patchi];
 243
 244             Info<< "Patch " << patchi
 245                 << " named " << nutPw.patch().name()
 246                 << " y+ : min: " << min(Yp) << " max: " << max(Yp)
 247                 << " average: " << average(Yp) << nl << endl;
 248         }
 249     }
 250
 251     if (!foundNutPatch)
 252     {
 253         Info<< "    no " << wallFunctionPatchField::typeName << " patches"
 254             << endl;
 255     }
 256 }
 257
 258
 259 int main(int argc, char *argv[])
 260 {
 261     timeSelector::addOptions();
 262
 263     #include "addRegionOption.H"
 264
 265     argList::validOptions.insert("compressible","");
 266     argList::validOptions.insert("twoPhase","");
 267
 268     #include "setRootCase.H"
 269     #include "createTime.H"
 270     instantList timeDirs = timeSelector::select0(runTime, args);
 271     #include "createNamedMesh.H"
 272
 273     bool compressible = args.optionFound("compressible");
 274
 275     // Check if two phase model was selected
 276     bool twoPhase = args.optionFound("twoPhase");
 277
 278     forAll(timeDirs, timeI)
 279     {
 280         runTime.setTime(timeDirs[timeI], timeI);
 281         Info<< "Time = " << runTime.timeName() << endl;
 282         fvMesh::readUpdateState state = mesh.readUpdate();
 283
 284         // Wall distance
 285         if (timeI == 0 || state != fvMesh::UNCHANGED)
 286         {
 287             Info<< "Calculating wall distance\n" << endl;
 288             wallDist y(mesh, true);
 289             Info<< "Writing wall distance to field "
 290                 << y.name() << nl << endl;
 291             y.write();
 292         }
 293
 294         volScalarField yPlus
 295         (
 296             IOobject
 297             (
 298                 "yPlus",
 299                 runTime.timeName(),
 300                 mesh,
 301                 IOobject::NO_READ,
 302                 IOobject::NO_WRITE
 303             ),
 304             mesh,
 305             dimensionedScalar("yPlus", dimless, 0.0)
 306         );
 307
 308         Info << "Reading field U\n" << endl;
 309         IOobject UHeader
 310         (
 311             "U",
 312             runTime.timeName(),
 313             mesh,
 314             IOobject::MUST_READ,
 315             IOobject::NO_WRITE
 316         );
 317
 318         if (UHeader.headerOk())
 319         {
 320             Info << "Reading field U\n" << endl;
 321             volVectorField U(UHeader, mesh);
 322
 323             if (compressible)
 324             {
 325                 calcCompressibleYPlus(mesh, runTime, U, yPlus);
 326             }
 327             else if (twoPhase)
 328             {
 329                 calcTwoPhaseYPlus(mesh, runTime, U, yPlus);
 330             }
 331             else
 332             {
 333                 calcIncompressibleYPlus(mesh, runTime, U, yPlus);
 334             }
 335         }
 336         else
 337         {
 338             Info<< "    no U field" << endl;
 339         }
 340
 341         Info<< "Writing yPlus to field " << yPlus.name() << nl << endl;
 342
 343         yPlus.write();
 344     }
 345
 346     Info<< "End\n" << endl;
 347
 348     return 0;
 349 }
 350
 351
 352 // ************************************************************************* //