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