src/turbulenceModels/RAS/incompressible/LRR/LRR.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 \*---------------------------------------------------------------------------*/
  26
  27 #include "LRR.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "wallFvPatch.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35 namespace incompressible
  36 {
  37 namespace RASModels
  38 {
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 defineTypeNameAndDebug(LRR, 0);
  43 addToRunTimeSelectionTable(RASModel, LRR, dictionary);
  44
  45 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  46
  47 // from components
  48 LRR::LRR
  49 (
  50     const volVectorField& U,
  51     const surfaceScalarField& phi,
  52     transportModel& lamTransportModel
  53 )
  54 :
  55     RASModel(typeName, U, phi, lamTransportModel),
  56
  57     Cmu_
  58     (
  59         dimensioned<scalar>::lookupOrAddToDict
  60         (
  61             "Cmu",
  62             coeffDict_,
  63             0.09
  64         )
  65     ),
  66     Clrr1_
  67     (
  68         dimensioned<scalar>::lookupOrAddToDict
  69         (
  70             "Clrr1",
  71             coeffDict_,
  72             1.8
  73         )
  74     ),
  75     Clrr2_
  76     (
  77         dimensioned<scalar>::lookupOrAddToDict
  78         (
  79             "Clrr2",
  80             coeffDict_,
  81             0.6
  82         )
  83     ),
  84     C1_
  85     (
  86         dimensioned<scalar>::lookupOrAddToDict
  87         (
  88             "C1",
  89             coeffDict_,
  90             1.44
  91         )
  92     ),
  93     C2_
  94     (
  95         dimensioned<scalar>::lookupOrAddToDict
  96         (
  97             "C2",
  98             coeffDict_,
  99             1.92
 100         )
 101     ),
 102     Cs_
 103     (
 104         dimensioned<scalar>::lookupOrAddToDict
 105         (
 106             "Cs",
 107             coeffDict_,
 108             0.25
 109         )
 110     ),
 111     Ceps_
 112     (
 113         dimensioned<scalar>::lookupOrAddToDict
 114         (
 115             "Ceps",
 116             coeffDict_,
 117             0.15
 118         )
 119     ),
 120     alphaEps_
 121     (
 122         dimensioned<scalar>::lookupOrAddToDict
 123         (
 124             "alphaEps",
 125             coeffDict_,
 126             0.76923
 127         )
 128     ),
 129     couplingFactor_
 130     (
 131         dimensioned<scalar>::lookupOrAddToDict
 132         (
 133             "couplingFactor",
 134             coeffDict_,
 135             0.0
 136         )
 137     ),
 138
 139     R_
 140     (
 141         IOobject
 142         (
 143             "R",
 144             runTime_.timeName(),
 145             mesh_,
 146             IOobject::MUST_READ,
 147             IOobject::AUTO_WRITE
 148         ),
 149         mesh_
 150     ),
 151
 152     k_
 153     (
 154         IOobject
 155         (
 156             "k",
 157             runTime_.timeName(),
 158             mesh_,
 159             IOobject::MUST_READ,
 160             IOobject::AUTO_WRITE
 161         ),
 162         mesh_
 163     ),
 164
 165     epsilon_
 166     (
 167         IOobject
 168         (
 169             "epsilon",
 170             runTime_.timeName(),
 171             mesh_,
 172             IOobject::MUST_READ,
 173             IOobject::AUTO_WRITE
 174         ),
 175         mesh_
 176     ),
 177
 178     nut_(Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_))
 179 {
 180 #   include "wallViscosityI.H"
 181
 182     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 183     {
 184         FatalErrorIn
 185         (
 186             "LRR::LRR"
 187             "(const volVectorField& U, const surfaceScalarField& phi,"
 188             "transportModel& lamTransportModel)"
 189         )   << "couplingFactor = " << couplingFactor_
 190             << " is not in range 0 - 1" << nl
 191             << exit(FatalError);
 192     }
 193
 194     printCoeffs();
 195 }
 196
 197
 198 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 199
 200 tmp<volSymmTensorField> LRR::devReff() const
 201 {
 202     return tmp<volSymmTensorField>
 203     (
 204         new volSymmTensorField
 205         (
 206             IOobject
 207             (
 208                 "devRhoReff",
 209                 runTime_.timeName(),
 210                 mesh_,
 211                 IOobject::NO_READ,
 212                 IOobject::NO_WRITE
 213             ),
 214             R_ - nu()*dev(twoSymm(fvc::grad(U_)))
 215         )
 216     );
 217 }
 218
 219
 220 tmp<fvVectorMatrix> LRR::divDevReff(volVectorField& U) const
 221 {
 222     if (couplingFactor_.value() > 0.0)
 223     {
 224         return
 225         (
 226             fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
 227           + fvc::laplacian
 228             (
 229                  (1.0 - couplingFactor_)*nut_,
 230                  U,
 231                  "laplacian(nuEff,U)"
 232             )
 233           - fvm::laplacian(nuEff(), U)
 234         );
 235     }
 236     else
 237     {
 238         return
 239         (
 240             fvc::div(R_)
 241           + fvc::laplacian(nut_, U, "laplacian(nuEff,U)")
 242           - fvm::laplacian(nuEff(), U)
 243         );
 244     }
 245 }
 246
 247
 248 bool LRR::read()
 249 {
 250     if (RASModel::read())
 251     {
 252         Cmu_.readIfPresent(coeffDict_);
 253         Clrr1_.readIfPresent(coeffDict_);
 254         Clrr2_.readIfPresent(coeffDict_);
 255         C1_.readIfPresent(coeffDict_);
 256         C2_.readIfPresent(coeffDict_);
 257         Cs_.readIfPresent(coeffDict_);
 258         Ceps_.readIfPresent(coeffDict_);
 259         alphaEps_.readIfPresent(coeffDict_);
 260
 261         couplingFactor_.readIfPresent(coeffDict_);
 262
 263         if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 264         {
 265             FatalErrorIn("LRR::read()")
 266                 << "couplingFactor = " << couplingFactor_
 267                 << " is not in range 0 - 1"
 268                 << exit(FatalError);
 269         }
 270
 271         return true;
 272     }
 273     else
 274     {
 275         return false;
 276     }
 277 }
 278
 279
 280 void LRR::correct()
 281 {
 282     transportModel_.correct();
 283
 284     if (!turbulence_)
 285     {
 286         return;
 287     }
 288
 289     RASModel::correct();
 290
 291     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
 292     volScalarField G = 0.5*mag(tr(P));
 293
 294 #   include "wallFunctionsI.H"
 295
 296     // Dissipation equation
 297     tmp<fvScalarMatrix> epsEqn
 298     (
 299         fvm::ddt(epsilon_)
 300       + fvm::div(phi_, epsilon_)
 301       - fvm::Sp(fvc::div(phi_), epsilon_)
 302     //- fvm::laplacian(Ceps*(K/epsilon_)*R, epsilon_)
 303       - fvm::laplacian(DepsilonEff(), epsilon_)
 304       ==
 305         C1_*G*epsilon_/k_
 306       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 307     );
 308
 309     epsEqn().relax();
 310
 311 #   include "wallDissipationI.H"
 312
 313     solve(epsEqn);
 314     bound(epsilon_, epsilon0_);
 315
 316
 317     // Reynolds stress equation
 318
 319     const fvPatchList& patches = mesh_.boundary();
 320
 321     forAll(patches, patchi)
 322     {
 323         const fvPatch& curPatch = patches[patchi];
 324
 325         if (typeid(curPatch) == typeid(wallFvPatch))
 326         {
 327             forAll(curPatch, facei)
 328             {
 329                 label faceCelli = curPatch.faceCells()[facei];
 330                 P[faceCelli]
 331                     *= min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 1.0);
 332             }
 333         }
 334     }
 335
 336
 337     tmp<fvSymmTensorMatrix> REqn
 338     (
 339         fvm::ddt(R_)
 340       + fvm::div(phi_, R_)
 341       - fvm::Sp(fvc::div(phi_), R_)
 342     //- fvm::laplacian(Cs*(k_/epsilon_)*R_, R_)
 343       - fvm::laplacian(DREff(), R_)
 344       + fvm::Sp(Clrr1_*epsilon_/k_, R_)
 345       ==
 346         P
 347       - (2.0/3.0*(1 - Clrr1_)*I)*epsilon_
 348       - Clrr2_*dev(P)
 349     );
 350
 351     REqn().relax();
 352     solve(REqn);
 353
 354     R_.max
 355     (
 356         dimensionedSymmTensor
 357         (
 358             "zero",
 359             R_.dimensions(),
 360             symmTensor
 361             (
 362                 k0_.value(), -GREAT, -GREAT,
 363                 k0_.value(), -GREAT,
 364                 k0_.value()
 365             )
 366         )
 367     );
 368
 369     k_ = 0.5*tr(R_);
 370     bound(k_, k0_);
 371
 372
 373     // Re-calculate viscosity
 374     nut_ = Cmu_*sqr(k_)/epsilon_;
 375
 376 #   include "wallViscosityI.H"
 377
 378
 379     // Correct wall shear stresses
 380
 381     forAll(patches, patchi)
 382     {
 383         const fvPatch& curPatch = patches[patchi];
 384
 385         if (typeid(curPatch) == typeid(wallFvPatch))
 386         {
 387             symmTensorField& Rw = R_.boundaryField()[patchi];
 388
 389             const scalarField& nutw = nut_.boundaryField()[patchi];
 390
 391             vectorField snGradU = U_.boundaryField()[patchi].snGrad();
 392
 393             const vectorField& faceAreas
 394                 = mesh_.Sf().boundaryField()[patchi];
 395
 396             const scalarField& magFaceAreas
 397                 = mesh_.magSf().boundaryField()[patchi];
 398
 399             forAll(curPatch, facei)
 400             {
 401                 // Calculate near-wall velocity gradient
 402                 tensor gradUw
 403                     = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
 404
 405                 // Calculate near-wall shear-stress tensor
 406                 tensor tauw = -nutw[facei]*2*symm(gradUw);
 407
 408                 // Reset the shear components of the stress tensor
 409                 Rw[facei].xy() = tauw.xy();
 410                 Rw[facei].xz() = tauw.xz();
 411                 Rw[facei].yz() = tauw.yz();
 412             }
 413         }
 414     }
 415 }
 416
 417
 418 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 419
 420 } // End namespace RASModels
 421 } // End namespace incompressible
 422 } // End namespace Foam
 423
 424 // ************************************************************************* //