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