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