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