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) 1991-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. To apply to
  32     compressible RAS cases, use the -compressible option.
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #include "fvCFD.H"
  37
  38 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
  39 #include "incompressible/RAS/RASModel/RASModel.H"
  40 #include "nutWallFunction/nutWallFunctionFvPatchScalarField.H"
  41
  42 #include "basicPsiThermo.H"
  43 #include "compressible/RAS/RASModel/RASModel.H"
  44 #include "mutWallFunction/mutWallFunctionFvPatchScalarField.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::nutWallFunctionFvPatchScalarField
  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::mutWallFunctionFvPatchScalarField
 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::validOptions.insert("compressible","");
 189
 190     #include "setRootCase.H"
 191     #include "createTime.H"
 192     instantList timeDirs = timeSelector::select0(runTime, args);
 193     #include "createNamedMesh.H"
 194
 195     bool compressible = args.optionFound("compressible");
 196
 197     forAll(timeDirs, timeI)
 198     {
 199         runTime.setTime(timeDirs[timeI], timeI);
 200         Info<< "Time = " << runTime.timeName() << endl;
 201         fvMesh::readUpdateState state = mesh.readUpdate();
 202
 203         // Wall distance
 204         if (timeI == 0 || state != fvMesh::UNCHANGED)
 205         {
 206             Info<< "Calculating wall distance\n" << endl;
 207             wallDist y(mesh, true);
 208             Info<< "Writing wall distance to field "
 209                 << y.name() << nl << endl;
 210             y.write();
 211         }
 212
 213         volScalarField yPlus
 214         (
 215             IOobject
 216             (
 217                 "yPlus",
 218                 runTime.timeName(),
 219                 mesh,
 220                 IOobject::NO_READ,
 221                 IOobject::NO_WRITE
 222             ),
 223             mesh,
 224             dimensionedScalar("yPlus", dimless, 0.0)
 225         );
 226
 227         IOobject UHeader
 228         (
 229             "U",
 230             runTime.timeName(),
 231             mesh,
 232             IOobject::MUST_READ,
 233             IOobject::NO_WRITE
 234         );
 235
 236         if (UHeader.headerOk())
 237         {
 238             Info << "Reading field U\n" << endl;
 239             volVectorField U(UHeader, mesh);
 240
 241             if (compressible)
 242             {
 243                 calcCompressibleYPlus(mesh, runTime, U, yPlus);
 244             }
 245             else
 246             {
 247                 calcIncompressibleYPlus(mesh, runTime, U, yPlus);
 248             }
 249         }
 250         else
 251         {
 252             Info<< "    no U field" << endl;
 253         }
 254
 255         Info<< "Writing yPlus to field " << yPlus.name() << nl << endl;
 256
 257         yPlus.write();
 258     }
 259
 260     Info<< "End\n" << endl;
 261
 262     return 0;
 263 }
 264
 265
 266 // ************************************************************************* //