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-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     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.optionFound("ybl"))
  80     {
  81         // If the boundary-layer thickness is provided use it
  82         ybl.value() = args.optionRead<scalar>("ybl");
  83     }
  84     else if (args.optionFound("Cbl"))
  85     {
  86         // Calculate boundary layer thickness as Cbl * mean distance to wall
  87         ybl.value() = gAverage(y) * args.optionRead<scalar>("Cbl");
  88     }
  89     else
  90     {
  91         FatalErrorIn(args.executable())
  92             << "Neither option 'ybl' or 'Cbl' have been provided to calculate"
  93                " the boundary-layer thickness"
  94             << exit(FatalError);
  95     }
  96
  97     Info<< "\nCreating boundary-layer for U of thickness "
  98         << ybl.value() << " m" << nl << endl;
  99
 100     // Modify velocity by applying a 1/7th power law boundary-layer
 101     // u/U0 = (y/ybl)^(1/7)
 102     // assumes U0 is the same as the current cell velocity
 103
 104     scalar yblv = ybl.value();
 105     forAll(U, celli)
 106     {
 107         if (y[celli] <= yblv)
 108         {
 109             U[celli] *= ::pow(y[celli]/yblv, (1.0/7.0));
 110         }
 111     }
 112
 113     Info<< "Writing U" << endl;
 114     U.write();
 115
 116     // Update/re-write phi
 117     phi = fvc::interpolate(U) & mesh.Sf();
 118     phi.write();
 119
 120     // Set turbulence constants
 121     dimensionedScalar kappa("kappa", dimless, 0.41);
 122     dimensionedScalar Cmu("Cmu", dimless, 0.09);
 123
 124     // Read and modify turbulence fields if present
 125
 126     IOobject epsilonHeader
 127     (
 128         "epsilon",
 129         runTime.timeName(),
 130         mesh,
 131         IOobject::MUST_READ
 132     );
 133
 134     IOobject kHeader
 135     (
 136         "k",
 137         runTime.timeName(),
 138         mesh,
 139         IOobject::MUST_READ
 140     );
 141
 142     IOobject nuTildaHeader
 143     (
 144         "nuTilda",
 145         runTime.timeName(),
 146         mesh,
 147         IOobject::MUST_READ
 148     );
 149
 150     // First calculate nut
 151     volScalarField nut
 152     (
 153         "nut",
 154         sqr(kappa*min(y, ybl))*::sqrt(2)*mag(dev(symm(fvc::grad(U))))
 155     );
 156
 157     if (args.optionFound("writenut"))
 158     {
 159         Info<< "Writing nut" << endl;
 160         nut.write();
 161     }
 162
 163
 164     // Read and modify turbulence fields if present
 165
 166     if (nuTildaHeader.headerOk())
 167     {
 168         Info<< "Reading field nuTilda\n" << endl;
 169         volScalarField nuTilda(nuTildaHeader, mesh);
 170         nuTilda = nut;
 171         nuTilda.correctBoundaryConditions();
 172
 173         Info<< "Writing nuTilda\n" << endl;
 174         nuTilda.write();
 175     }
 176
 177     if (kHeader.headerOk() && epsilonHeader.headerOk())
 178     {
 179         Info<< "Reading field k\n" << endl;
 180         volScalarField k(kHeader, mesh);
 181
 182         Info<< "Reading field epsilon\n" << endl;
 183         volScalarField epsilon(epsilonHeader, mesh);
 184
 185         scalar ck0 = ::pow(Cmu.value(), 0.25)*kappa.value();
 186         k = sqr(nut/(ck0*min(y, ybl)));
 187         k.correctBoundaryConditions();
 188
 189         scalar ce0 = ::pow(Cmu.value(), 0.75)/kappa.value();
 190         epsilon = ce0*k*sqrt(k)/min(y, ybl);
 191         epsilon.correctBoundaryConditions();
 192
 193         Info<< "Writing k\n" << endl;
 194         k.write();
 195
 196         Info<< "Writing epsilon\n" << endl;
 197         epsilon.write();
 198     }
 199
 200     Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
 201         << "  ClockTime = " << runTime.elapsedClockTime() << " s"
 202         << nl << endl;
 203
 204     Info<< "End\n" << endl;
 205
 206     return 0;
 207 }
 208
 209
 210 // ************************************************************************* //