applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.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     applyBoundaryLayer
  26
  27 Description
  28     Apply a simplified boundary-layer model to the velocity and
  29     turbulence fields based on the 1/7th power-law.
  30
  31     The uniform boundary-layer thickness is either provided via the -ybl option
  32     or calculated as the average of the distance to the wall scaled with
  33     the thickness coefficient supplied via the option -Cbl.  If both options
  34     are provided -ybl is used.
  35
  36 \*---------------------------------------------------------------------------*/
  37
  38 #include "fvCFD.H"
  39 #include "singlePhaseTransportModel.H"
  40 #include "RASModel.H"
  41 #include "wallDist.H"
  42
  43 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  44
  45 // turbulence constants - file-scope
  46 static const scalar Cmu(0.09);
  47 static const scalar kappa(0.41);
  48
  49
  50 int main(int argc, char *argv[])
  51 {
  52     argList::addNote
  53     (
  54         "apply a simplified boundary-layer model to the velocity and\n"
  55         "turbulence fields based on the 1/7th power-law."
  56     );
  57
  58     argList::addOption
  59     (
  60         "ybl",
  61         "scalar",
  62         "specify the boundary-layer thickness"
  63     );
  64     argList::addOption
  65     (
  66         "Cbl", "scalar",
  67         "boundary-layer thickness as Cbl * mean distance to wall"
  68     );
  69     argList::addBoolOption
  70     (
  71         "writenut",
  72         "write nut field"
  73     );
  74
  75     #include "setRootCase.H"
  76     #include "createTime.H"
  77     #include "createMesh.H"
  78     #include "createFields.H"
  79
  80     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  81
  82     Info<< "Time = " << runTime.timeName() << nl << endl;
  83
  84     // Modify velocity by applying a 1/7th power law boundary-layer
  85     // u/U0 = (y/ybl)^(1/7)
  86     // assumes U0 is the same as the current cell velocity
  87
  88     scalar yblv = ybl.value();
  89     forAll(U, celli)
  90     {
  91         if (y[celli] <= yblv)
  92         {
  93             U[celli] *= ::pow(y[celli]/yblv, (1.0/7.0));
  94         }
  95     }
  96
  97     Info<< "Writing U\n" << endl;
  98     U.write();
  99
 100     // Update/re-write phi
 101     phi = fvc::interpolate(U) & mesh.Sf();
 102     phi.write();
 103
 104     // Calculate nut
 105     tmp<volScalarField> tnut = turbulence->nut();
 106     volScalarField& nut = tnut();
 107     volScalarField S(mag(dev(symm(fvc::grad(U)))));
 108     nut = sqr(kappa*min(y, ybl))*::sqrt(2)*S;
 109
 110     if (args.optionFound("writenut"))
 111     {
 112         Info<< "Writing nut" << endl;
 113         nut.write();
 114     }
 115
 116     // Create G field - used by RAS wall functions
 117     volScalarField G("RASModel::G", nut*2*sqr(S));
 118
 119
 120     //--- Read and modify turbulence fields
 121
 122     // Turbulence k
 123     tmp<volScalarField> tk = turbulence->k();
 124     volScalarField& k = tk();
 125     scalar ck0 = pow025(Cmu)*kappa;
 126     k = sqr(nut/(ck0*min(y, ybl)));
 127     k.correctBoundaryConditions();
 128
 129     Info<< "Writing k\n" << endl;
 130     k.write();
 131
 132
 133     // Turbulence epsilon
 134     tmp<volScalarField> tepsilon = turbulence->epsilon();
 135     volScalarField& epsilon = tepsilon();
 136     scalar ce0 = ::pow(Cmu, 0.75)/kappa;
 137     epsilon = ce0*k*sqrt(k)/min(y, ybl);
 138     epsilon.correctBoundaryConditions();
 139
 140     Info<< "Writing epsilon\n" << endl;
 141     epsilon.write();
 142
 143
 144     // Turbulence omega
 145     IOobject omegaHeader
 146     (
 147         "omega",
 148         runTime.timeName(),
 149         mesh,
 150         IOobject::MUST_READ,
 151         IOobject::NO_WRITE,
 152         false
 153     );
 154
 155     if (omegaHeader.headerOk())
 156     {
 157         volScalarField omega(omegaHeader, mesh);
 158         omega = epsilon/(Cmu*k);
 159         omega.correctBoundaryConditions();
 160
 161         Info<< "Writing omega\n" << endl;
 162         omega.write();
 163     }
 164
 165
 166     // Turbulence nuTilda
 167     IOobject nuTildaHeader
 168     (
 169         "nuTilda",
 170         runTime.timeName(),
 171         mesh,
 172         IOobject::MUST_READ,
 173         IOobject::NO_WRITE,
 174         false
 175     );
 176
 177     if (nuTildaHeader.headerOk())
 178     {
 179         volScalarField nuTilda(nuTildaHeader, mesh);
 180         nuTilda = nut;
 181         nuTilda.correctBoundaryConditions();
 182
 183         Info<< "Writing nuTilda\n" << endl;
 184         nuTilda.write();
 185     }
 186
 187
 188     Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 189         << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 190         << nl << endl;
 191
 192     Info<< "End\n" << endl;
 193
 194     return 0;
 195 }
 196
 197
 198 // ************************************************************************* //