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