src/turbulenceModels/incompressible/RAS/LRR/LRR.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "LRR.H"
  27 #include "addToRunTimeSelectionTable.H"
  28
  29 #include "backwardsCompatibilityWallFunctions.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 LRR::LRR
  48 (
  49     const volVectorField& U,
  50     const surfaceScalarField& phi,
  51     transportModel& lamTransportModel
  52 )
  53 :
  54     RASModel(typeName, U, phi, lamTransportModel),
  55
  56     Cmu_
  57     (
  58         dimensionedScalar::lookupOrAddToDict
  59         (
  60             "Cmu",
  61             coeffDict_,
  62             0.09
  63         )
  64     ),
  65     Clrr1_
  66     (
  67         dimensionedScalar::lookupOrAddToDict
  68         (
  69             "Clrr1",
  70             coeffDict_,
  71             1.8
  72         )
  73     ),
  74     Clrr2_
  75     (
  76         dimensionedScalar::lookupOrAddToDict
  77         (
  78             "Clrr2",
  79             coeffDict_,
  80             0.6
  81         )
  82     ),
  83     C1_
  84     (
  85         dimensionedScalar::lookupOrAddToDict
  86         (
  87             "C1",
  88             coeffDict_,
  89             1.44
  90         )
  91     ),
  92     C2_
  93     (
  94         dimensionedScalar::lookupOrAddToDict
  95         (
  96             "C2",
  97             coeffDict_,
  98             1.92
  99         )
 100     ),
 101     Cs_
 102     (
 103         dimensionedScalar::lookupOrAddToDict
 104         (
 105             "Cs",
 106             coeffDict_,
 107             0.25
 108         )
 109     ),
 110     Ceps_
 111     (
 112         dimensionedScalar::lookupOrAddToDict
 113         (
 114             "Ceps",
 115             coeffDict_,
 116             0.15
 117         )
 118     ),
 119     sigmaEps_
 120     (
 121         dimensionedScalar::lookupOrAddToDict
 122         (
 123             "sigmaEps",
 124             coeffDict_,
 125             1.3
 126         )
 127     ),
 128     couplingFactor_
 129     (
 130         dimensionedScalar::lookupOrAddToDict
 131         (
 132             "couplingFactor",
 133             coeffDict_,
 134             0.0
 135         )
 136     ),
 137
 138     R_
 139     (
 140         IOobject
 141         (
 142             "R",
 143             runTime_.timeName(),
 144             U_.db(),
 145             IOobject::NO_READ,
 146             IOobject::AUTO_WRITE
 147         ),
 148         autoCreateR("R", mesh_)
 149     ),
 150     k_
 151     (
 152         IOobject
 153         (
 154             "k",
 155             runTime_.timeName(),
 156             U_.db(),
 157             IOobject::NO_READ,
 158             IOobject::AUTO_WRITE
 159         ),
 160         autoCreateK("k", mesh_)
 161     ),
 162     epsilon_
 163     (
 164         IOobject
 165         (
 166             "epsilon",
 167             runTime_.timeName(),
 168             U_.db(),
 169             IOobject::NO_READ,
 170             IOobject::AUTO_WRITE
 171         ),
 172         autoCreateEpsilon("epsilon", mesh_)
 173     ),
 174     nut_
 175     (
 176         IOobject
 177         (
 178             "nut",
 179             runTime_.timeName(),
 180             U_.db(),
 181             IOobject::NO_READ,
 182             IOobject::AUTO_WRITE
 183         ),
 184         autoCreateNut("nut", mesh_)
 185     )
 186 {
 187     nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
 188     nut_ = min(nut_, nuRatio()*nu());
 189     nut_.correctBoundaryConditions();
 190
 191     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 192     {
 193         FatalErrorIn
 194         (
 195             "LRR::LRR"
 196             "(const volVectorField& U, const surfaceScalarField& phi,"
 197             "transportModel& lamTransportModel)"
 198         )   << "couplingFactor = " << couplingFactor_
 199             << " is not in range 0 - 1" << nl
 200             << exit(FatalError);
 201     }
 202
 203     printCoeffs();
 204 }
 205
 206
 207 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 208
 209 tmp<volSymmTensorField> LRR::devReff() const
 210 {
 211     return tmp<volSymmTensorField>
 212     (
 213         new volSymmTensorField
 214         (
 215             IOobject
 216             (
 217                 "devRhoReff",
 218                 runTime_.timeName(),
 219                 U_.db(),
 220                 IOobject::NO_READ,
 221                 IOobject::NO_WRITE
 222             ),
 223             R_ - nu()*dev(twoSymm(fvc::grad(U_)))
 224         )
 225     );
 226 }
 227
 228
 229 tmp<fvVectorMatrix> LRR::divDevReff(volVectorField& U) const
 230 {
 231     if (couplingFactor_.value() > 0.0)
 232     {
 233         return
 234         (
 235             fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
 236           + fvc::laplacian
 237             (
 238                  (1.0 - couplingFactor_)*nut_,
 239                  U,
 240                  "laplacian(nuEff,U)"
 241             )
 242           - fvm::laplacian(nuEff(), U)
 243         );
 244     }
 245     else
 246     {
 247         return
 248         (
 249             fvc::div(R_)
 250           + fvc::laplacian(nut_, U, "laplacian(nuEff,U)")
 251           - fvm::laplacian(nuEff(), U)
 252         );
 253     }
 254 }
 255
 256
 257 bool LRR::read()
 258 {
 259     if (RASModel::read())
 260     {
 261         Cmu_.readIfPresent(coeffDict());
 262         Clrr1_.readIfPresent(coeffDict());
 263         Clrr2_.readIfPresent(coeffDict());
 264         C1_.readIfPresent(coeffDict());
 265         C2_.readIfPresent(coeffDict());
 266         Cs_.readIfPresent(coeffDict());
 267         Ceps_.readIfPresent(coeffDict());
 268         sigmaEps_.readIfPresent(coeffDict());
 269
 270         couplingFactor_.readIfPresent(coeffDict());
 271
 272         if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 273         {
 274             FatalErrorIn("LRR::read()")
 275                 << "couplingFactor = " << couplingFactor_
 276                 << " is not in range 0 - 1"
 277                 << exit(FatalError);
 278         }
 279
 280         return true;
 281     }
 282     else
 283     {
 284         return false;
 285     }
 286 }
 287
 288
 289 void LRR::correct()
 290 {
 291     // Bound in case of topological change
 292     // HJ, 22/Aug/2007
 293     if (mesh_.changing())
 294     {
 295         R_.correctBoundaryConditions();
 296         bound(epsilon_, epsilon0_);
 297     }
 298
 299     RASModel::correct();
 300
 301     if (!turbulence_)
 302     {
 303         return;
 304     }
 305
 306     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
 307     volScalarField G("RASModel::G", 0.5*mag(tr(P)));
 308
 309     // Update epsilon and G at the wall
 310     epsilon_.boundaryField().updateCoeffs();
 311
 312     // Dissipation equation
 313     tmp<fvScalarMatrix> epsEqn
 314     (
 315         fvm::ddt(epsilon_)
 316       + fvm::div(phi_, epsilon_)
 317       + fvm::SuSp(-fvc::div(phi_), epsilon_)
 318     //- fvm::laplacian(Ceps*(K/epsilon_)*R, epsilon_)
 319       - fvm::laplacian(DepsilonEff(), epsilon_)
 320       ==
 321         C1_*G*epsilon_/k_
 322       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 323     );
 324
 325     epsEqn().relax();
 326
 327     // No longer needed: matrix completes at the point of solution
 328     // HJ, 17/Apr/2012
 329 //     epsEqn().completeAssembly();
 330
 331     solve(epsEqn);
 332     bound(epsilon_, epsilon0_);
 333
 334
 335     // Reynolds stress equation
 336
 337     const fvPatchList& patches = mesh_.boundary();
 338
 339     forAll(patches, patchi)
 340     {
 341         const fvPatch& curPatch = patches[patchi];
 342
 343         if (curPatch.isWall())
 344         {
 345             forAll(curPatch, facei)
 346             {
 347                 label faceCelli = curPatch.faceCells()[facei];
 348
 349                 P[faceCelli] *=
 350                     min
 351                     (
 352                         G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL),
 353                         1.0
 354                     );
 355             }
 356         }
 357     }
 358
 359
 360     tmp<fvSymmTensorMatrix> REqn
 361     (
 362         fvm::ddt(R_)
 363       + fvm::div(phi_, R_)
 364       + fvm::SuSp(-fvc::div(phi_), R_)
 365     //- fvm::laplacian(Cs*(k_/epsilon_)*R_, R_)
 366       - fvm::laplacian(DREff(), R_)
 367       + fvm::Sp(Clrr1_*epsilon_/k_, R_)
 368       ==
 369         P
 370       - (2.0/3.0*(1 - Clrr1_)*I)*epsilon_
 371       - Clrr2_*dev(P)
 372     );
 373
 374     REqn().relax();
 375     solve(REqn);
 376
 377     R_.max
 378     (
 379         dimensionedSymmTensor
 380         (
 381             "zero",
 382             R_.dimensions(),
 383             symmTensor
 384             (
 385                 k0_.value(), -GREAT, -GREAT,
 386                 k0_.value(), -GREAT,
 387                 k0_.value()
 388             )
 389         )
 390     );
 391
 392     k_ = 0.5*tr(R_);
 393     bound(k_, k0_);
 394
 395
 396     // Re-calculate viscosity
 397     nut_ = Cmu_*sqr(k_)/epsilon_;
 398     nut_ = min(nut_, nuRatio()*nu());
 399     nut_.correctBoundaryConditions();
 400
 401
 402     // Correct wall shear stresses
 403
 404     forAll(patches, patchi)
 405     {
 406         const fvPatch& curPatch = patches[patchi];
 407
 408         if (curPatch.isWall())
 409         {
 410             symmTensorField& Rw = R_.boundaryField()[patchi];
 411
 412             const scalarField& nutw = nut_.boundaryField()[patchi];
 413
 414             vectorField snGradU = U_.boundaryField()[patchi].snGrad();
 415
 416             const vectorField& faceAreas
 417                 = mesh_.Sf().boundaryField()[patchi];
 418
 419             const scalarField& magFaceAreas
 420                 = mesh_.magSf().boundaryField()[patchi];
 421
 422             forAll(curPatch, facei)
 423             {
 424                 // Calculate near-wall velocity gradient
 425                 tensor gradUw
 426                     = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
 427
 428                 // Calculate near-wall shear-stress tensor
 429                 tensor tauw = -nutw[facei]*2*symm(gradUw);
 430
 431                 // Reset the shear components of the stress tensor
 432                 Rw[facei].xy() = tauw.xy();
 433                 Rw[facei].xz() = tauw.xz();
 434                 Rw[facei].yz() = tauw.yz();
 435             }
 436         }
 437     }
 438 }
 439
 440
 441 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 442
 443 } // End namespace RASModels
 444 } // End namespace incompressible
 445 } // End namespace Foam
 446
 447 // ************************************************************************* //