src/turbulenceModels/RAS/compressible/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 #include "wallDist.H"
  31 #include "wallDistReflection.H"
  32
  33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  34
  35 namespace Foam
  36 {
  37 namespace compressible
  38 {
  39 namespace RASModels
  40 {
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 defineTypeNameAndDebug(LaunderGibsonRSTM, 0);
  45 addToRunTimeSelectionTable(RASModel, LaunderGibsonRSTM, dictionary);
  46
  47 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  48
  49 // from components
  50 LaunderGibsonRSTM::LaunderGibsonRSTM
  51 (
  52     const volScalarField& rho,
  53     const volVectorField& U,
  54     const surfaceScalarField& phi,
  55     basicThermo& thermophysicalModel
  56 )
  57 :
  58     RASModel(typeName, rho, U, phi, thermophysicalModel),
  59
  60     Cmu_
  61     (
  62         dimensioned<scalar>::lookupOrAddToDict
  63         (
  64             "Cmu",
  65             coeffDict_,
  66             0.09
  67         )
  68     ),
  69     Clg1_
  70     (
  71         dimensioned<scalar>::lookupOrAddToDict
  72         (
  73             "Clg1",
  74             coeffDict_,
  75             1.8
  76         )
  77     ),
  78     Clg2_
  79     (
  80         dimensioned<scalar>::lookupOrAddToDict
  81         (
  82             "Clg2",
  83             coeffDict_,
  84             0.6
  85         )
  86     ),
  87     C1_
  88     (
  89         dimensioned<scalar>::lookupOrAddToDict
  90         (
  91             "C1",
  92             coeffDict_,
  93             1.44
  94         )
  95     ),
  96     C2_
  97     (
  98         dimensioned<scalar>::lookupOrAddToDict
  99         (
 100             "C2",
 101             coeffDict_,
 102             1.92
 103         )
 104     ),
 105     Cs_
 106     (
 107         dimensioned<scalar>::lookupOrAddToDict
 108         (
 109             "Cs",
 110             coeffDict_,
 111             0.25
 112         )
 113     ),
 114     Ceps_
 115     (
 116         dimensioned<scalar>::lookupOrAddToDict
 117         (
 118             "Ceps",
 119             coeffDict_,
 120             0.15
 121         )
 122     ),
 123     C1Ref_
 124     (
 125         dimensioned<scalar>::lookupOrAddToDict
 126         (
 127             "C1Ref",
 128             coeffDict_,
 129             0.5
 130         )
 131     ),
 132     C2Ref_
 133     (
 134         dimensioned<scalar>::lookupOrAddToDict
 135         (
 136             "C2Ref",
 137             coeffDict_,
 138             0.3
 139         )
 140     ),
 141     couplingFactor_
 142     (
 143         dimensioned<scalar>::lookupOrAddToDict
 144         (
 145             "couplingFactor",
 146             coeffDict_,
 147             0.0
 148         )
 149     ),
 150     alphaR_
 151     (
 152         dimensioned<scalar>::lookupOrAddToDict
 153         (
 154             "alphaR",
 155             coeffDict_,
 156             1.22
 157         )
 158     ),
 159     alphaEps_
 160     (
 161         dimensioned<scalar>::lookupOrAddToDict
 162         (
 163             "alphaEps",
 164             coeffDict_,
 165             0.76923
 166         )
 167     ),
 168     alphah_
 169     (
 170         dimensioned<scalar>::lookupOrAddToDict
 171         (
 172             "alphah",
 173             coeffDict_,
 174             1.0
 175         )
 176     ),
 177
 178     y_(mesh_),
 179
 180     R_
 181     (
 182         IOobject
 183         (
 184             "R",
 185             runTime_.timeName(),
 186             mesh_,
 187             IOobject::MUST_READ,
 188             IOobject::AUTO_WRITE
 189         ),
 190         mesh_
 191     ),
 192
 193     k_
 194     (
 195         IOobject
 196         (
 197             "k",
 198             runTime_.timeName(),
 199             mesh_,
 200             IOobject::MUST_READ,
 201             IOobject::AUTO_WRITE
 202         ),
 203         mesh_
 204     ),
 205
 206     epsilon_
 207     (
 208         IOobject
 209         (
 210             "epsilon",
 211             runTime_.timeName(),
 212             mesh_,
 213             IOobject::MUST_READ,
 214             IOobject::AUTO_WRITE
 215         ),
 216         mesh_
 217     ),
 218
 219     mut_
 220     (
 221         IOobject
 222         (
 223             "mut",
 224             runTime_.timeName(),
 225             mesh_,
 226             IOobject::NO_READ,
 227             IOobject::NO_WRITE
 228         ),
 229         Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_)
 230     )
 231 {
 232 #   include "wallViscosityI.H"
 233
 234     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 235     {
 236         FatalErrorIn
 237         (
 238             "LaunderGibsonRSTM::LaunderGibsonRSTM"
 239             "(const volScalarField&, const volVectorField&"
 240             ", const surfaceScalarField&, basicThermo&)"
 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::devRhoReff() 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             rho_*R_ - mu()*dev(twoSymm(fvc::grad(U_)))
 267         )
 268     );
 269 }
 270
 271
 272 tmp<fvVectorMatrix> LaunderGibsonRSTM::divDevRhoReff(volVectorField& U) const
 273 {
 274     if (couplingFactor_.value() > 0.0)
 275     {
 276         return
 277         (
 278             fvc::div(rho_*R_ + couplingFactor_*mut_*fvc::grad(U))
 279           + fvc::laplacian((1.0 - couplingFactor_)*mut_, U)
 280           - fvm::laplacian(muEff(), U)
 281           - fvc::div(mu()*dev2(fvc::grad(U)().T()))
 282         );
 283     }
 284     else
 285     {
 286         return
 287         (
 288             fvc::div(rho_*R_)
 289           + fvc::laplacian(mut_, U)
 290           - fvm::laplacian(muEff(), U)
 291           - fvc::div(mu()*dev2(fvc::grad(U)().T()))
 292         );
 293     }
 294 }
 295
 296
 297 bool LaunderGibsonRSTM::read()
 298 {
 299     if (RASModel::read())
 300     {
 301         Cmu_.readIfPresent(coeffDict_);
 302         Clg1_.readIfPresent(coeffDict_);
 303         Clg2_.readIfPresent(coeffDict_);
 304         C1_.readIfPresent(coeffDict_);
 305         C2_.readIfPresent(coeffDict_);
 306         Cs_.readIfPresent(coeffDict_);
 307         Ceps_.readIfPresent(coeffDict_);
 308         C1Ref_.readIfPresent(coeffDict_);
 309         C2Ref_.readIfPresent(coeffDict_);
 310         alphaR_.readIfPresent(coeffDict_);
 311         alphaEps_.readIfPresent(coeffDict_);
 312         alphah_.readIfPresent(coeffDict_);
 313
 314         couplingFactor_.readIfPresent(coeffDict_);
 315
 316         if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 317         {
 318             FatalErrorIn("LaunderGibsonRSTM::read()")
 319                 << "couplingFactor = " << couplingFactor_
 320                 << " is not in range 0 - 1" << nl
 321                 << exit(FatalError);
 322         }
 323
 324         return true;
 325     }
 326     else
 327     {
 328         return false;
 329     }
 330 }
 331
 332
 333 void LaunderGibsonRSTM::correct()
 334 {
 335     if (!turbulence_)
 336     {
 337         // Re-calculate viscosity
 338         mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
 339         return;
 340     }
 341
 342     RASModel::correct();
 343
 344     if (mesh_.changing())
 345     {
 346         y_.correct();
 347     }
 348
 349     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
 350     volScalarField G = 0.5*mag(tr(P));
 351
 352 #   include "wallFunctionsI.H"
 353
 354     // Dissipation equation
 355     tmp<fvScalarMatrix> epsEqn
 356     (
 357         fvm::ddt(rho_, epsilon_)
 358       + fvm::div(phi_, epsilon_)
 359     //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
 360       - fvm::laplacian(DepsilonEff(), epsilon_)
 361      ==
 362         C1_*rho_*G*epsilon_/k_
 363       - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
 364     );
 365
 366     epsEqn().relax();
 367
 368 #   include "wallDissipationI.H"
 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 (typeid(curPatch) == typeid(wallFvPatch))
 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), 100.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(rho_, R_)
 398       + fvm::div(phi_, R_)
 399     //- fvm::laplacian(Cs*rho_*(k_/epsilon_)*R_, R_)
 400       - fvm::laplacian(DREff(), R_)
 401       + fvm::Sp(Clg1_*rho_*epsilon_/k_, R_)
 402      ==
 403         rho_*P
 404       + (2.0/3.0*(Clg1_ - 1)*I)*rho_*epsilon_
 405       - Clg2_*rho_*dev(P)
 406
 407         // wall reflection terms
 408       + symm
 409         (
 410             I*((y_.n() & reflect) & y_.n())
 411           - 1.5*(y_.n()*(reflect & y_.n())
 412           + (y_.n() & reflect)*y_.n())
 413         )*pow(Cmu_, 0.75)*rho_*pow(k_, 1.5)/(kappa_*y_*epsilon_)
 414     );
 415
 416     REqn().relax();
 417     solve(REqn);
 418
 419     R_.max
 420     (
 421         dimensionedSymmTensor
 422         (
 423             "zero",
 424             R_.dimensions(),
 425             symmTensor
 426             (
 427                 k0_.value(), -GREAT, -GREAT,
 428                 k0_.value(), -GREAT,
 429                 k0_.value()
 430             )
 431         )
 432     );
 433
 434     k_ == 0.5*tr(R_);
 435     bound(k_, k0_);
 436
 437
 438     // Re-calculate turbulent viscosity
 439     mut_ = Cmu_*rho_*sqr(k_)/epsilon_;
 440
 441
 442 #   include "wallViscosityI.H"
 443
 444
 445     // Correct wall shear stresses
 446
 447     forAll(patches, patchi)
 448     {
 449         const fvPatch& curPatch = patches[patchi];
 450
 451         if (typeid(curPatch) == typeid(wallFvPatch))
 452         {
 453             symmTensorField& Rw = R_.boundaryField()[patchi];
 454
 455             const scalarField& mutw = mut_.boundaryField()[patchi];
 456             const scalarField& rhow = rho_.boundaryField()[patchi];
 457
 458             vectorField snGradU = U_.boundaryField()[patchi].snGrad();
 459
 460             const vectorField& faceAreas
 461                 = mesh_.Sf().boundaryField()[patchi];
 462
 463             const scalarField& magFaceAreas
 464                 = mesh_.magSf().boundaryField()[patchi];
 465
 466             forAll(curPatch, facei)
 467             {
 468                 // Calculate near-wall velocity gradient
 469                 tensor gradUw
 470                     = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
 471
 472                 // Calculate near-wall shear-stress tensor
 473                 tensor tauw = -(mutw[facei]/rhow[facei])*2*dev(symm(gradUw));
 474
 475                 // Reset the shear components of the stress tensor
 476                 Rw[facei].xy() = tauw.xy();
 477                 Rw[facei].xz() = tauw.xz();
 478                 Rw[facei].yz() = tauw.yz();
 479             }
 480         }
 481     }
 482 }
 483
 484
 485 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 486
 487 } // End namespace RASModels
 488 } // End namespace compressible
 489 } // End namespace Foam
 490
 491 // ************************************************************************* //