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