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 held by original author
   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
  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& 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     sigmaEps_
 121     (
 122         dimensioned<scalar>::lookupOrAddToDict
 123         (
 124             "sigmaEps",
 125             coeffDict_,
 126             1.3
 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::NO_READ,
 147             IOobject::AUTO_WRITE
 148         ),
 149         autoCreateR("R", mesh_)
 150     ),
 151     k_
 152     (
 153         IOobject
 154         (
 155             "k",
 156             runTime_.timeName(),
 157             mesh_,
 158             IOobject::NO_READ,
 159             IOobject::AUTO_WRITE
 160         ),
 161         autoCreateK("k", mesh_)
 162     ),
 163     epsilon_
 164     (
 165         IOobject
 166         (
 167             "epsilon",
 168             runTime_.timeName(),
 169             mesh_,
 170             IOobject::NO_READ,
 171             IOobject::AUTO_WRITE
 172         ),
 173         autoCreateEpsilon("epsilon", mesh_)
 174     ),
 175     nut_
 176     (
 177         IOobject
 178         (
 179             "nut",
 180             runTime_.timeName(),
 181             mesh_,
 182             IOobject::NO_READ,
 183             IOobject::AUTO_WRITE
 184         ),
 185         autoCreateNut("nut", mesh_)
 186     )
 187 {
 188     nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
 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                 mesh_,
 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     RASModel::correct();
 292
 293     if (!turbulence_)
 294     {
 295         return;
 296     }
 297
 298     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
 299     volScalarField G("RASModel::G", 0.5*mag(tr(P)));
 300
 301     // Update espsilon and G at the wall
 302     epsilon_.boundaryField().updateCoeffs();
 303
 304     // Dissipation equation
 305     tmp<fvScalarMatrix> epsEqn
 306     (
 307         fvm::ddt(epsilon_)
 308       + fvm::div(phi_, epsilon_)
 309       + fvm::SuSp(-fvc::div(phi_), epsilon_)
 310     //- fvm::laplacian(Ceps*(K/epsilon_)*R, epsilon_)
 311       - fvm::laplacian(DepsilonEff(), epsilon_)
 312       ==
 313         C1_*G*epsilon_/k_
 314       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 315     );
 316
 317     epsEqn().relax();
 318
 319     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 320
 321     solve(epsEqn);
 322     bound(epsilon_, epsilon0_);
 323
 324
 325     // Reynolds stress equation
 326
 327     const fvPatchList& patches = mesh_.boundary();
 328
 329     forAll(patches, patchi)
 330     {
 331         const fvPatch& curPatch = patches[patchi];
 332
 333         if (curPatch.isWall())
 334         {
 335             forAll(curPatch, facei)
 336             {
 337                 label faceCelli = curPatch.faceCells()[facei];
 338                 P[faceCelli]
 339                     *= min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 1.0);
 340             }
 341         }
 342     }
 343
 344
 345     tmp<fvSymmTensorMatrix> REqn
 346     (
 347         fvm::ddt(R_)
 348       + fvm::div(phi_, R_)
 349       + fvm::SuSp(-fvc::div(phi_), R_)
 350     //- fvm::laplacian(Cs*(k_/epsilon_)*R_, R_)
 351       - fvm::laplacian(DREff(), R_)
 352       + fvm::Sp(Clrr1_*epsilon_/k_, R_)
 353       ==
 354         P
 355       - (2.0/3.0*(1 - Clrr1_)*I)*epsilon_
 356       - Clrr2_*dev(P)
 357     );
 358
 359     REqn().relax();
 360     solve(REqn);
 361
 362     R_.max
 363     (
 364         dimensionedSymmTensor
 365         (
 366             "zero",
 367             R_.dimensions(),
 368             symmTensor
 369             (
 370                 k0_.value(), -GREAT, -GREAT,
 371                 k0_.value(), -GREAT,
 372                 k0_.value()
 373             )
 374         )
 375     );
 376
 377     k_ = 0.5*tr(R_);
 378     bound(k_, k0_);
 379
 380
 381     // Re-calculate viscosity
 382     nut_ = Cmu_*sqr(k_)/epsilon_;
 383     nut_.correctBoundaryConditions();
 384
 385
 386     // Correct wall shear stresses
 387
 388     forAll(patches, patchi)
 389     {
 390         const fvPatch& curPatch = patches[patchi];
 391
 392         if (curPatch.isWall())
 393         {
 394             symmTensorField& Rw = R_.boundaryField()[patchi];
 395
 396             const scalarField& nutw = nut_.boundaryField()[patchi];
 397
 398             vectorField snGradU = U_.boundaryField()[patchi].snGrad();
 399
 400             const vectorField& faceAreas
 401                 = mesh_.Sf().boundaryField()[patchi];
 402
 403             const scalarField& magFaceAreas
 404                 = mesh_.magSf().boundaryField()[patchi];
 405
 406             forAll(curPatch, facei)
 407             {
 408                 // Calculate near-wall velocity gradient
 409                 tensor gradUw
 410                     = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
 411
 412                 // Calculate near-wall shear-stress tensor
 413                 tensor tauw = -nutw[facei]*2*symm(gradUw);
 414
 415                 // Reset the shear components of the stress tensor
 416                 Rw[facei].xy() = tauw.xy();
 417                 Rw[facei].xz() = tauw.xz();
 418                 Rw[facei].yz() = tauw.yz();
 419             }
 420         }
 421     }
 422 }
 423
 424
 425 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 426
 427 } // End namespace RASModels
 428 } // End namespace incompressible
 429 } // End namespace Foam
 430
 431 // ************************************************************************* //