src/turbulenceModels/compressible/RAS/LRR/LRR.C

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