src/turbulenceModels/compressible/RAS/LaunderGibsonRSTM/LaunderGibsonRSTM.C

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