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