src/turbulenceModels/RAS/incompressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 "LaunderGibsonRSTM.H"
  28 #include "addToRunTimeSelectionTable.H"
  29 #include "wallFvPatch.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35 namespace incompressible
  36 {
  37 namespace RASModels
  38 {
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 defineTypeNameAndDebug(LaunderGibsonRSTM, 0);
  43 addToRunTimeSelectionTable(RASModel, LaunderGibsonRSTM, dictionary);
  44
  45 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  46
  47 // from components
  48 LaunderGibsonRSTM::LaunderGibsonRSTM
  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     Clg1_
  67     (
  68         dimensioned<scalar>::lookupOrAddToDict
  69         (
  70             "Clg1",
  71             coeffDict_,
  72             1.8
  73         )
  74     ),
  75     Clg2_
  76     (
  77         dimensioned<scalar>::lookupOrAddToDict
  78         (
  79             "Clg2",
  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     alphaR_
 121     (
 122         dimensioned<scalar>::lookupOrAddToDict
 123         (
 124             "alphaR",
 125             coeffDict_,
 126             1.22
 127         )
 128     ),
 129     alphaEps_
 130     (
 131         dimensioned<scalar>::lookupOrAddToDict
 132         (
 133             "alphaEps",
 134             coeffDict_,
 135             0.76923
 136         )
 137     ),
 138     C1Ref_
 139     (
 140         dimensioned<scalar>::lookupOrAddToDict
 141         (
 142             "C1Ref",
 143             coeffDict_,
 144             0.5
 145         )
 146     ),
 147     C2Ref_
 148     (
 149         dimensioned<scalar>::lookupOrAddToDict
 150         (
 151             "C2Ref",
 152             coeffDict_,
 153             0.3
 154         )
 155     ),
 156     couplingFactor_
 157     (
 158         dimensioned<scalar>::lookupOrAddToDict
 159         (
 160             "couplingFactor",
 161             coeffDict_,
 162             0.0
 163         )
 164     ),
 165
 166     yr_(mesh_),
 167
 168     R_
 169     (
 170         IOobject
 171         (
 172             "R",
 173             runTime_.timeName(),
 174             mesh_,
 175             IOobject::MUST_READ,
 176             IOobject::AUTO_WRITE
 177         ),
 178         mesh_
 179     ),
 180
 181     k_
 182     (
 183         IOobject
 184         (
 185             "k",
 186             runTime_.timeName(),
 187             mesh_,
 188             IOobject::MUST_READ,
 189             IOobject::AUTO_WRITE
 190         ),
 191         mesh_
 192     ),
 193
 194     epsilon_
 195     (
 196         IOobject
 197         (
 198             "epsilon",
 199             runTime_.timeName(),
 200             mesh_,
 201             IOobject::MUST_READ,
 202             IOobject::AUTO_WRITE
 203         ),
 204         mesh_
 205     ),
 206
 207     nut_(Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_))
 208 {
 209 #   include "wallViscosityI.H"
 210
 211     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 212     {
 213         FatalErrorIn
 214         (
 215             "LaunderGibsonRSTM::LaunderGibsonRSTM"
 216             "(const volVectorField& U, const surfaceScalarField& phi,"
 217             "transportModel& lamTransportModel)"
 218         )   << "couplingFactor = " << couplingFactor_
 219             << " is not in range 0 - 1" << nl
 220             << exit(FatalError);
 221     }
 222
 223     printCoeffs();
 224 }
 225
 226
 227 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 228
 229 tmp<volSymmTensorField> LaunderGibsonRSTM::devReff() const
 230 {
 231     return tmp<volSymmTensorField>
 232     (
 233         new volSymmTensorField
 234         (
 235             IOobject
 236             (
 237                 "devRhoReff",
 238                 runTime_.timeName(),
 239                 mesh_,
 240                 IOobject::NO_READ,
 241                 IOobject::NO_WRITE
 242             ),
 243             R_ - nu()*dev(twoSymm(fvc::grad(U_)))
 244         )
 245     );
 246 }
 247
 248
 249 tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevReff(volVectorField& U) const
 250 {
 251     if (couplingFactor_.value() > 0.0)
 252     {
 253         return
 254         (
 255             fvc::div(R_ + couplingFactor_*nut_*fvc::grad(U), "div(R)")
 256           + fvc::laplacian((1.0-couplingFactor_)*nut_, U, "laplacian(nuEff,U)")
 257           - fvm::laplacian(nuEff(), U)
 258         );
 259     }
 260     else
 261     {
 262         return
 263         (
 264             fvc::div(R_)
 265           + fvc::laplacian(nut_, U, "laplacian(nuEff,U)")
 266           - fvm::laplacian(nuEff(), U)
 267         );
 268     }
 269 }
 270
 271
 272 bool LaunderGibsonRSTM::read()
 273 {
 274     if (RASModel::read())
 275     {
 276         Cmu_.readIfPresent(coeffDict_);
 277         Clg1_.readIfPresent(coeffDict_);
 278         Clg2_.readIfPresent(coeffDict_);
 279         C1_.readIfPresent(coeffDict_);
 280         C2_.readIfPresent(coeffDict_);
 281         Cs_.readIfPresent(coeffDict_);
 282         Ceps_.readIfPresent(coeffDict_);
 283         alphaR_.readIfPresent(coeffDict_);
 284         alphaEps_.readIfPresent(coeffDict_);
 285         C1Ref_.readIfPresent(coeffDict_);
 286         C2Ref_.readIfPresent(coeffDict_);
 287
 288         couplingFactor_.readIfPresent(coeffDict_);
 289
 290         if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 291         {
 292             FatalErrorIn("LaunderGibsonRSTM::read()")
 293                 << "couplingFactor = " << couplingFactor_
 294                 << " is not in range 0 - 1"
 295                 << exit(FatalError);
 296         }
 297
 298         return true;
 299     }
 300     else
 301     {
 302         return false;
 303     }
 304 }
 305
 306
 307 void LaunderGibsonRSTM::correct()
 308 {
 309     transportModel_.correct();
 310
 311     if (!turbulence_)
 312     {
 313         return;
 314     }
 315
 316     RASModel::correct();
 317
 318     if (mesh_.changing())
 319     {
 320         yr_.correct();
 321     }
 322
 323     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
 324     volScalarField G = 0.5*tr(P);
 325
 326 #   include "wallFunctionsI.H"
 327
 328     // Dissipation equation
 329     tmp<fvScalarMatrix> epsEqn
 330     (
 331         fvm::ddt(epsilon_)
 332       + fvm::div(phi_, epsilon_)
 333     //- fvm::laplacian(Ceps*(k_/epsilon_)*R_, epsilon_)
 334       - fvm::laplacian(DepsilonEff(), epsilon_)
 335      ==
 336         C1_*G*epsilon_/k_
 337       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 338     );
 339
 340     epsEqn().relax();
 341
 342 #   include "wallDissipationI.H"
 343
 344     solve(epsEqn);
 345     bound(epsilon_, epsilon0_);
 346
 347
 348     // Reynolds stress equation
 349
 350     const fvPatchList& patches = mesh_.boundary();
 351
 352     forAll(patches, patchi)
 353     {
 354         const fvPatch& curPatch = patches[patchi];
 355
 356         if (typeid(curPatch) == typeid(wallFvPatch))
 357         {
 358             forAll(curPatch, facei)
 359             {
 360                 label faceCelli = curPatch.faceCells()[facei];
 361                 P[faceCelli] *=
 362                     min(G[faceCelli]/(0.5*tr(P[faceCelli]) + SMALL), 1.0);
 363             }
 364         }
 365     }
 366
 367     volSymmTensorField reflect = C1Ref_*epsilon_/k_*R_ - C2Ref_*Clg2_*dev(P);
 368
 369     tmp<fvSymmTensorMatrix> REqn
 370     (
 371         fvm::ddt(R_)
 372       + fvm::div(phi_, R_)
 373     //- fvm::laplacian(Cs*(k_/epsilon_)*R_, R_)
 374       - fvm::laplacian(DREff(), R_)
 375       + fvm::Sp(Clg1_*epsilon_/k_, R_)
 376       ==
 377         P
 378       + (2.0/3.0*(Clg1_ - 1)*I)*epsilon_
 379       - Clg2_*dev(P)
 380
 381         // wall reflection terms
 382       + symm
 383         (
 384             I*((yr_.n() & reflect) & yr_.n())
 385           - 1.5*(yr_.n()*(reflect & yr_.n())
 386           + (yr_.n() & reflect)*yr_.n())
 387         )*pow(Cmu_, 0.75)*pow(k_, 1.5)/(kappa_*yr_*epsilon_)
 388     );
 389
 390     REqn().relax();
 391     solve(REqn);
 392
 393     R_.max
 394     (
 395         dimensionedSymmTensor
 396         (
 397             "zero",
 398             R_.dimensions(),
 399             symmTensor
 400             (
 401                 k0_.value(), -GREAT, -GREAT,
 402                              k0_.value(), -GREAT,
 403                                           k0_.value()
 404             )
 405         )
 406     );
 407
 408     k_ == 0.5*tr(R_);
 409     bound(k_, k0_);
 410
 411
 412     // Re-calculate turbulent viscosity
 413     nut_ = Cmu_*sqr(k_)/epsilon_;
 414
 415
 416 #   include "wallViscosityI.H"
 417
 418
 419     // Correct wall shear stresses
 420
 421     forAll(patches, patchi)
 422     {
 423         const fvPatch& curPatch = patches[patchi];
 424
 425         if (typeid(curPatch) == typeid(wallFvPatch))
 426         {
 427             symmTensorField& Rw = R_.boundaryField()[patchi];
 428
 429             const scalarField& nutw = nut_.boundaryField()[patchi];
 430
 431             vectorField snGradU = U_.boundaryField()[patchi].snGrad();
 432
 433             const vectorField& faceAreas
 434                 = mesh_.Sf().boundaryField()[patchi];
 435
 436             const scalarField& magFaceAreas
 437                 = mesh_.magSf().boundaryField()[patchi];
 438
 439             forAll(curPatch, facei)
 440             {
 441                 // Calculate near-wall velocity gradient
 442                 tensor gradUw
 443                     = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
 444
 445                 // Calculate near-wall shear-stress tensor
 446                 tensor tauw = -nutw[facei]*2*symm(gradUw);
 447
 448                 // Reset the shear components of the stress tensor
 449                 Rw[facei].xy() = tauw.xy();
 450                 Rw[facei].xz() = tauw.xz();
 451                 Rw[facei].yz() = tauw.yz();
 452             }
 453         }
 454     }
 455 }
 456
 457
 458 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 459
 460 } // End namespace RASModels
 461 } // End namespace incompressible
 462 } // End namespace Foam
 463
 464 // ************************************************************************* //