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 (C) 2004-2010 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, 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.
  32     Use the -compressible option for compressible RAS cases.
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "fvCFD.H"
  37
  38 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
  39 #include "incompressible/RAS/RASModel/RASModel.H"
  40 #include "nutkWallFunction/nutkWallFunctionFvPatchScalarField.H"
  41
  42 #include "basicPsiThermo.H"
  43 #include "compressible/RAS/RASModel/RASModel.H"
  44 #include "mutkWallFunction/mutkWallFunctionFvPatchScalarField.H"
  45
  46 #include "wallDist.H"
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49
  50 void calcIncompressibleYPlus
  51 (
  52     const fvMesh& mesh,
  53     const Time& runTime,
  54     const volVectorField& U,
  55     volScalarField& yPlus
  56 )
  57 {
  58     typedef incompressible::RASModels::nutkWallFunctionFvPatchScalarField
  59         wallFunctionPatchField;
  60
  61     #include "createPhi.H"
  62
  63     singlePhaseTransportModel laminarTransport(U, phi);
  64
  65     autoPtr<incompressible::RASModel> RASModel
  66     (
  67         incompressible::RASModel::New(U, phi, laminarTransport)
  68     );
  69
  70     const volScalarField::GeometricBoundaryField nutPatches =
  71         RASModel->nut()().boundaryField();
  72
  73     bool foundNutPatch = false;
  74     forAll(nutPatches, patchi)
  75     {
  76         if (isA<wallFunctionPatchField>(nutPatches[patchi]))
  77         {
  78             foundNutPatch = true;
  79
  80             const wallFunctionPatchField& nutPw =
  81                 dynamic_cast<const wallFunctionPatchField&>
  82                     (nutPatches[patchi]);
  83
  84             yPlus.boundaryField()[patchi] = nutPw.yPlus();
  85             const scalarField& Yp = yPlus.boundaryField()[patchi];
  86
  87             Info<< "Patch " << patchi
  88                 << " named " << nutPw.patch().name()
  89                 << " y+ : min: " << min(Yp) << " max: " << max(Yp)
  90                 << " average: " << average(Yp) << nl << endl;
  91         }
  92     }
  93
  94     if (!foundNutPatch)
  95     {
  96         Info<< "    no " << wallFunctionPatchField::typeName << " patches"
  97             << endl;
  98     }
  99 }
 100
 101
 102 void calcCompressibleYPlus
 103 (
 104     const fvMesh& mesh,
 105     const Time& runTime,
 106     const volVectorField& U,
 107     volScalarField& yPlus
 108 )
 109 {
 110     typedef compressible::RASModels::mutkWallFunctionFvPatchScalarField
 111         wallFunctionPatchField;
 112
 113     IOobject rhoHeader
 114     (
 115         "rho",
 116         runTime.timeName(),
 117         mesh,
 118         IOobject::MUST_READ,
 119         IOobject::NO_WRITE
 120     );
 121
 122     if (!rhoHeader.headerOk())
 123     {
 124         Info<< "    no rho field" << endl;
 125         return;
 126     }
 127
 128     Info<< "Reading field rho\n" << endl;
 129     volScalarField rho(rhoHeader, mesh);
 130
 131     #include "compressibleCreatePhi.H"
 132
 133     autoPtr<basicPsiThermo> pThermo
 134     (
 135         basicPsiThermo::New(mesh)
 136     );
 137     basicPsiThermo& thermo = pThermo();
 138
 139     autoPtr<compressible::RASModel> RASModel
 140     (
 141         compressible::RASModel::New
 142         (
 143             rho,
 144             U,
 145             phi,
 146             thermo
 147         )
 148     );
 149
 150     const volScalarField::GeometricBoundaryField mutPatches =
 151         RASModel->mut()().boundaryField();
 152
 153     bool foundMutPatch = false;
 154     forAll(mutPatches, patchi)
 155     {
 156         if (isA<wallFunctionPatchField>(mutPatches[patchi]))
 157         {
 158             foundMutPatch = true;
 159
 160             const wallFunctionPatchField& mutPw =
 161                 dynamic_cast<const wallFunctionPatchField&>
 162                     (mutPatches[patchi]);
 163
 164             yPlus.boundaryField()[patchi] = mutPw.yPlus();
 165             const scalarField& Yp = yPlus.boundaryField()[patchi];
 166
 167             Info<< "Patch " << patchi
 168                 << " named " << mutPw.patch().name()
 169                 << " y+ : min: " << min(Yp) << " max: " << max(Yp)
 170                 << " average: " << average(Yp) << nl << endl;
 171         }
 172     }
 173
 174     if (!foundMutPatch)
 175     {
 176         Info<< "    no " << wallFunctionPatchField::typeName << " patches"
 177             << endl;
 178     }
 179 }
 180
 181
 182 int main(int argc, char *argv[])
 183 {
 184     timeSelector::addOptions();
 185
 186     #include "addRegionOption.H"
 187
 188     argList::addBoolOption
 189     (
 190         "compressible",
 191         "calculate compressible y+"
 192     );
 193
 194     #include "setRootCase.H"
 195     #include "createTime.H"
 196     instantList timeDirs = timeSelector::select0(runTime, args);
 197     #include "createNamedMesh.H"
 198
 199     const bool compressible = args.optionFound("compressible");
 200
 201     forAll(timeDirs, timeI)
 202     {
 203         runTime.setTime(timeDirs[timeI], timeI);
 204         Info<< "Time = " << runTime.timeName() << endl;
 205         fvMesh::readUpdateState state = mesh.readUpdate();
 206
 207         // Wall distance
 208         if (timeI == 0 || state != fvMesh::UNCHANGED)
 209         {
 210             Info<< "Calculating wall distance\n" << endl;
 211             wallDist y(mesh, true);
 212             Info<< "Writing wall distance to field " << y.name() << nl << endl;
 213             y.write();
 214         }
 215
 216         volScalarField yPlus
 217         (
 218             IOobject
 219             (
 220                 "yPlus",
 221                 runTime.timeName(),
 222                 mesh,
 223                 IOobject::NO_READ,
 224                 IOobject::NO_WRITE
 225             ),
 226             mesh,
 227             dimensionedScalar("yPlus", dimless, 0.0)
 228         );
 229
 230         IOobject UHeader
 231         (
 232             "U",
 233             runTime.timeName(),
 234             mesh,
 235             IOobject::MUST_READ,
 236             IOobject::NO_WRITE
 237         );
 238
 239         if (UHeader.headerOk())
 240         {
 241             Info<< "Reading field U\n" << endl;
 242             volVectorField U(UHeader, mesh);
 243
 244             if (compressible)
 245             {
 246                 calcCompressibleYPlus(mesh, runTime, U, yPlus);
 247             }
 248             else
 249             {
 250                 calcIncompressibleYPlus(mesh, runTime, U, yPlus);
 251             }
 252         }
 253         else
 254         {
 255             Info<< "    no U field" << endl;
 256         }
 257
 258         Info<< "Writing yPlus to field " << yPlus.name() << nl << endl;
 259
 260         yPlus.write();
 261     }
 262
 263     Info<< "End\n" << endl;
 264
 265     return 0;
 266 }
 267
 268
 269 // ************************************************************************* //