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 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 "LRR.H"
  28 #include "addToRunTimeSelectionTable.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 )
  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     Clrr1_
  68     (
  69         dimensioned<scalar>::lookupOrAddToDict
  70         (
  71             "Clrr1",
  72             coeffDict_,
  73             1.8
  74         )
  75     ),
  76     Clrr2_
  77     (
  78         dimensioned<scalar>::lookupOrAddToDict
  79         (
  80             "Clrr2",
  81             coeffDict_,
  82             0.6
  83         )
  84     ),
  85     C1_
  86     (
  87         dimensioned<scalar>::lookupOrAddToDict
  88         (
  89             "C1",
  90             coeffDict_,
  91             1.44
  92         )
  93     ),
  94     C2_
  95     (
  96         dimensioned<scalar>::lookupOrAddToDict
  97         (
  98             "C2",
  99             coeffDict_,
 100             1.92
 101         )
 102     ),
 103     Cs_
 104     (
 105         dimensioned<scalar>::lookupOrAddToDict
 106         (
 107             "Cs",
 108             coeffDict_,
 109             0.25
 110         )
 111     ),
 112     Ceps_
 113     (
 114         dimensioned<scalar>::lookupOrAddToDict
 115         (
 116             "Ceps",
 117             coeffDict_,
 118             0.15
 119         )
 120     ),
 121     couplingFactor_
 122     (
 123         dimensioned<scalar>::lookupOrAddToDict
 124         (
 125             "couplingFactor",
 126             coeffDict_,
 127             0.0
 128         )
 129     ),
 130     sigmaR_
 131     (
 132         dimensioned<scalar>::lookupOrAddToDict
 133         (
 134             "sigmaR",
 135             coeffDict_,
 136             0.81967
 137         )
 138     ),
 139     sigmaEps_
 140     (
 141         dimensioned<scalar>::lookupOrAddToDict
 142         (
 143             "sigmaEps",
 144             coeffDict_,
 145             1.3
 146         )
 147     ),
 148     Prt_
 149     (
 150         dimensioned<scalar>::lookupOrAddToDict
 151         (
 152             "Prt",
 153             coeffDict_,
 154             1.0
 155         )
 156     ),
 157
 158     R_
 159     (
 160         IOobject
 161         (
 162             "R",
 163             runTime_.timeName(),
 164             mesh_,
 165             IOobject::MUST_READ,
 166             IOobject::AUTO_WRITE
 167         ),
 168         autoCreateR("R", mesh_)
 169     ),
 170     k_
 171     (
 172         IOobject
 173         (
 174             "k",
 175             runTime_.timeName(),
 176             mesh_,
 177             IOobject::NO_READ,
 178             IOobject::AUTO_WRITE
 179         ),
 180         autoCreateK("k", mesh_)
 181     ),
 182     epsilon_
 183     (
 184         IOobject
 185         (
 186             "epsilon",
 187             runTime_.timeName(),
 188             mesh_,
 189             IOobject::NO_READ,
 190             IOobject::AUTO_WRITE
 191         ),
 192         autoCreateEpsilon("epsilon", mesh_)
 193     ),
 194     mut_
 195     (
 196         IOobject
 197         (
 198             "mut",
 199             runTime_.timeName(),
 200             mesh_,
 201             IOobject::NO_READ,
 202             IOobject::AUTO_WRITE
 203         ),
 204         autoCreateMut("mut", mesh_)
 205     ),
 206     alphat_
 207     (
 208         IOobject
 209         (
 210             "alphat",
 211             runTime_.timeName(),
 212             mesh_,
 213             IOobject::NO_READ,
 214             IOobject::AUTO_WRITE
 215         ),
 216         autoCreateAlphat("alphat", mesh_)
 217     )
 218 {
 219     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 220     {
 221         FatalErrorIn
 222         (
 223             "LRR::LRR"
 224             "( const volScalarField&, const volVectorField&"
 225             ", const surfaceScalarField&, incompressibleTransportModel&)"
 226         )   << "couplingFactor = " << couplingFactor_
 227             << " is not in range 0 - 1" << nl
 228             << exit(FatalError);
 229     }
 230
 231     mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
 232     mut_.correctBoundaryConditions();
 233
 234     alphat_ = mut_/Prt_;
 235     alphat_.correctBoundaryConditions();
 236
 237     printCoeffs();
 238 }
 239
 240
 241 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 242
 243 tmp<volSymmTensorField> LRR::devRhoReff() const
 244 {
 245     return tmp<volSymmTensorField>
 246     (
 247         new volSymmTensorField
 248         (
 249             IOobject
 250             (
 251                 "devRhoReff",
 252                 runTime_.timeName(),
 253                 mesh_,
 254                 IOobject::NO_READ,
 255                 IOobject::NO_WRITE
 256             ),
 257             rho_*R_ - mu()*dev(twoSymm(fvc::grad(U_)))
 258         )
 259     );
 260 }
 261
 262
 263 tmp<fvVectorMatrix> LRR::divDevRhoReff(volVectorField& U) const
 264 {
 265     if (couplingFactor_.value() > 0.0)
 266     {
 267         return
 268         (
 269             fvc::div(rho_*R_ + couplingFactor_*mut_*fvc::grad(U))
 270           + fvc::laplacian((1.0 - couplingFactor_)*mut_, U)
 271           - fvm::laplacian(muEff(), U)
 272           - fvc::div(mu()*dev2(fvc::grad(U)().T()))
 273         );
 274     }
 275     else
 276     {
 277         return
 278         (
 279             fvc::div(rho_*R_)
 280           + fvc::laplacian(mut_, U)
 281           - fvm::laplacian(muEff(), U)
 282           - fvc::div(mu()*dev2(fvc::grad(U)().T()))
 283         );
 284     }
 285 }
 286
 287
 288 bool LRR::read()
 289 {
 290     if (RASModel::read())
 291     {
 292         Cmu_.readIfPresent(coeffDict());
 293         Clrr1_.readIfPresent(coeffDict());
 294         Clrr2_.readIfPresent(coeffDict());
 295         C1_.readIfPresent(coeffDict());
 296         C2_.readIfPresent(coeffDict());
 297         Cs_.readIfPresent(coeffDict());
 298         Ceps_.readIfPresent(coeffDict());
 299         sigmaR_.readIfPresent(coeffDict());
 300         sigmaEps_.readIfPresent(coeffDict());
 301         Prt_.readIfPresent(coeffDict());
 302         couplingFactor_.readIfPresent(coeffDict());
 303
 304         if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
 305         {
 306             FatalErrorIn("LRR::read()")
 307                 << "couplingFactor = " << couplingFactor_
 308                 << " is not in range 0 - 1" << nl
 309                 << exit(FatalError);
 310         }
 311
 312         return true;
 313     }
 314     else
 315     {
 316         return false;
 317     }
 318 }
 319
 320
 321 void LRR::correct()
 322 {
 323     // Bound in case of topological change
 324     // HJ, 22/Aug/2007
 325     if (mesh_.changing())
 326     {
 327         bound(k_, k0_);
 328         bound(epsilon_, epsilon0_);
 329     }
 330
 331     if (!turbulence_)
 332     {
 333         // Re-calculate viscosity
 334         mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
 335         mut_.correctBoundaryConditions();
 336
 337         // Re-calculate thermal diffusivity
 338         alphat_ = mut_/Prt_;
 339         alphat_.correctBoundaryConditions();
 340
 341         return;
 342     }
 343
 344     RASModel::correct();
 345
 346     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
 347     volScalarField G("RASModel::G", 0.5*mag(tr(P)));
 348
 349     // Update epsilon and G at the wall
 350     epsilon_.boundaryField().updateCoeffs();
 351
 352     // Dissipation equation
 353     tmp<fvScalarMatrix> epsEqn
 354     (
 355         fvm::ddt(rho_, epsilon_)
 356       + fvm::div(phi_, epsilon_)
 357     //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
 358       - fvm::laplacian(DepsilonEff(), epsilon_)
 359      ==
 360         C1_*rho_*G*epsilon_/k_
 361       - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
 362     );
 363
 364     epsEqn().relax();
 365
 366     epsEqn().boundaryManipulate(epsilon_.boundaryField());
 367
 368     solve(epsEqn);
 369     bound(epsilon_, epsilon0_);
 370
 371
 372     // Reynolds stress equation
 373
 374     const fvPatchList& patches = mesh_.boundary();
 375
 376     forAll(patches, patchi)
 377     {
 378         const fvPatch& curPatch = patches[patchi];
 379
 380         if (curPatch.isWall())
 381         {
 382             forAll(curPatch, facei)
 383             {
 384                 label faceCelli = curPatch.faceCells()[facei];
 385                 P[faceCelli]
 386                     *= min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 100.0);
 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                 k0_.value(), -GREAT, -GREAT,
 417                 k0_.value(), -GREAT,
 418                 k0_.value()
 419             )
 420         )
 421     );
 422
 423     k_ = 0.5*tr(R_);
 424     bound(k_, k0_);
 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 (curPatch.isWall())
 443         {
 444             symmTensorField& Rw = R_.boundaryField()[patchi];
 445
 446             const scalarField& rhow = rho_.boundaryField()[patchi];
 447             const scalarField& mutw = mut_.boundaryField()[patchi];
 448
 449             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 // ************************************************************************* //