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