src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 "NonlinearKEShih.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "wallFvPatch.H"
  29
  30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  31
  32 namespace Foam
  33 {
  34 namespace incompressible
  35 {
  36 namespace RASModels
  37 {
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 defineTypeNameAndDebug(NonlinearKEShih, 0);
  42 addToRunTimeSelectionTable(RASModel, NonlinearKEShih, dictionary);
  43
  44 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  45
  46 NonlinearKEShih::NonlinearKEShih
  47 (
  48     const volVectorField& U,
  49     const surfaceScalarField& phi,
  50     transportModel& transport,
  51     const word& turbulenceModelName,
  52     const word& modelName
  53 )
  54 :
  55     RASModel(modelName, U, phi, transport, turbulenceModelName),
  56
  57     C1_
  58     (
  59         dimensioned<scalar>::lookupOrAddToDict
  60         (
  61             "C1",
  62             coeffDict_,
  63             1.44
  64         )
  65     ),
  66     C2_
  67     (
  68         dimensioned<scalar>::lookupOrAddToDict
  69         (
  70             "C2",
  71             coeffDict_,
  72             1.92
  73         )
  74     ),
  75     sigmak_
  76     (
  77         dimensioned<scalar>::lookupOrAddToDict
  78         (
  79             "sigmak",
  80             coeffDict_,
  81             1.0
  82         )
  83     ),
  84     sigmaEps_
  85     (
  86         dimensioned<scalar>::lookupOrAddToDict
  87         (
  88             "sigmaEps",
  89             coeffDict_,
  90             1.3
  91         )
  92     ),
  93     A1_
  94     (
  95         dimensioned<scalar>::lookupOrAddToDict
  96         (
  97             "A1",
  98             coeffDict_,
  99             1.25
 100         )
 101     ),
 102     A2_
 103     (
 104         dimensioned<scalar>::lookupOrAddToDict
 105         (
 106             "A2",
 107             coeffDict_,
 108             1000.0
 109         )
 110     ),
 111     Ctau1_
 112     (
 113         dimensioned<scalar>::lookupOrAddToDict
 114         (
 115             "Ctau1",
 116             coeffDict_,
 117             -4.0
 118         )
 119     ),
 120     Ctau2_
 121     (
 122         dimensioned<scalar>::lookupOrAddToDict
 123         (
 124             "Ctau2",
 125             coeffDict_,
 126             13.0
 127         )
 128     ),
 129     Ctau3_
 130     (
 131         dimensioned<scalar>::lookupOrAddToDict
 132         (
 133             "Ctau3",
 134             coeffDict_,
 135             -2.0
 136         )
 137     ),
 138     alphaKsi_
 139     (
 140         dimensioned<scalar>::lookupOrAddToDict
 141         (
 142             "alphaKsi",
 143             coeffDict_,
 144             0.9
 145         )
 146     ),
 147
 148     kappa_
 149     (
 150         dimensioned<scalar>::lookupOrAddToDict
 151         (
 152             "kappa_",
 153             coeffDict_,
 154             0.41
 155         )
 156     ),
 157     E_
 158     (
 159         dimensioned<scalar>::lookupOrAddToDict
 160         (
 161             "E",
 162             coeffDict_,
 163             9.8
 164         )
 165     ),
 166
 167     k_
 168     (
 169         IOobject
 170         (
 171             "k",
 172             runTime_.timeName(),
 173             mesh_,
 174             IOobject::MUST_READ,
 175             IOobject::AUTO_WRITE
 176         ),
 177         mesh_
 178     ),
 179
 180     epsilon_
 181     (
 182         IOobject
 183         (
 184             "epsilon",
 185             runTime_.timeName(),
 186             mesh_,
 187             IOobject::MUST_READ,
 188             IOobject::AUTO_WRITE
 189         ),
 190         mesh_
 191     ),
 192
 193     gradU_(fvc::grad(U)),
 194     eta_
 195     (
 196         k_/bound(epsilon_, epsilonMin_)
 197        *sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())))
 198     ),
 199     ksi_
 200     (
 201         k_/epsilon_
 202        *sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())))
 203     ),
 204     Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
 205     fEta_(A2_ + pow(eta_, 3.0)),
 206
 207     nut_("nut", Cmu_*sqr(k_)/epsilon_),
 208
 209     nonlinearStress_
 210     (
 211         "nonlinearStress",
 212         symm
 213         (
 214             pow3(k_)/sqr(epsilon_)
 215            *(
 216                 Ctau1_/fEta_
 217                *(
 218                     (gradU_ & gradU_)
 219                   + (gradU_ & gradU_)().T()
 220                 )
 221               + Ctau2_/fEta_*(gradU_ & gradU_.T())
 222               + Ctau3_/fEta_*(gradU_.T() & gradU_)
 223             )
 224         )
 225     )
 226 {
 227     bound(k_, kMin_);
 228     // already bounded: bound(epsilon_, epsilonMin_);
 229
 230     #include "wallNonlinearViscosityI.H"
 231
 232     printCoeffs();
 233 }
 234
 235
 236 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 237
 238 tmp<volSymmTensorField> NonlinearKEShih::R() const
 239 {
 240     return tmp<volSymmTensorField>
 241     (
 242         new volSymmTensorField
 243         (
 244             IOobject
 245             (
 246                 "R",
 247                 runTime_.timeName(),
 248                 mesh_,
 249                 IOobject::NO_READ,
 250                 IOobject::NO_WRITE
 251             ),
 252             ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)) + nonlinearStress_,
 253             k_.boundaryField().types()
 254         )
 255     );
 256 }
 257
 258
 259 tmp<volSymmTensorField> NonlinearKEShih::devReff() const
 260 {
 261     return tmp<volSymmTensorField>
 262     (
 263         new volSymmTensorField
 264         (
 265             IOobject
 266             (
 267                 "devRhoReff",
 268                 runTime_.timeName(),
 269                 mesh_,
 270                 IOobject::NO_READ,
 271                 IOobject::NO_WRITE
 272             ),
 273            -nuEff()*dev(twoSymm(fvc::grad(U_))) + nonlinearStress_
 274         )
 275     );
 276 }
 277
 278
 279 tmp<fvVectorMatrix> NonlinearKEShih::divDevReff(volVectorField& U) const
 280 {
 281     return
 282     (
 283         fvc::div(nonlinearStress_)
 284       - fvm::laplacian(nuEff(), U)
 285       - fvc::div(nuEff()*dev(T(fvc::grad(U))))
 286     );
 287 }
 288
 289
 290 bool NonlinearKEShih::read()
 291 {
 292     if (RASModel::read())
 293     {
 294         C1_.readIfPresent(coeffDict());
 295         C2_.readIfPresent(coeffDict());
 296         sigmak_.readIfPresent(coeffDict());
 297         sigmaEps_.readIfPresent(coeffDict());
 298         A1_.readIfPresent(coeffDict());
 299         A2_.readIfPresent(coeffDict());
 300         Ctau1_.readIfPresent(coeffDict());
 301         Ctau2_.readIfPresent(coeffDict());
 302         Ctau3_.readIfPresent(coeffDict());
 303         alphaKsi_.readIfPresent(coeffDict());
 304
 305         kappa_.readIfPresent(coeffDict());
 306         E_.readIfPresent(coeffDict());
 307
 308         return true;
 309     }
 310     else
 311     {
 312         return false;
 313     }
 314 }
 315
 316
 317 void NonlinearKEShih::correct()
 318 {
 319     RASModel::correct();
 320
 321     if (!turbulence_)
 322     {
 323         return;
 324     }
 325
 326     gradU_ = fvc::grad(U_);
 327
 328     // generation term
 329     tmp<volScalarField> S2 = symm(gradU_) && gradU_;
 330
 331     volScalarField G
 332     (
 333         "RASModel::G",
 334         Cmu_*sqr(k_)/epsilon_*S2
 335       - (nonlinearStress_ && gradU_)
 336     );
 337
 338     #include "nonLinearWallFunctionsI.H"
 339
 340     // Dissipation equation
 341     tmp<fvScalarMatrix> epsEqn
 342     (
 343         fvm::ddt(epsilon_)
 344       + fvm::div(phi_, epsilon_)
 345       - fvm::laplacian(DepsilonEff(), epsilon_)
 346       ==
 347         C1_*G*epsilon_/k_
 348       - fvm::Sp(C2_*epsilon_/k_, epsilon_)
 349     );
 350
 351     epsEqn().relax();
 352
 353     #include "wallDissipationI.H"
 354
 355     solve(epsEqn);
 356     bound(epsilon_, epsilonMin_);
 357
 358
 359     // Turbulent kinetic energy equation
 360
 361     tmp<fvScalarMatrix> kEqn
 362     (
 363         fvm::ddt(k_)
 364       + fvm::div(phi_, k_)
 365       - fvm::laplacian(DkEff(), k_)
 366       ==
 367         G
 368       - fvm::Sp(epsilon_/k_, k_)
 369     );
 370
 371     kEqn().relax();
 372     solve(kEqn);
 373     bound(k_, kMin_);
 374
 375
 376     // Re-calculate viscosity
 377
 378     eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ + gradU_.T())));
 379     ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU_ - gradU_.T())));
 380     Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
 381     fEta_ = A2_ + pow(eta_, 3.0);
 382
 383     nut_ = Cmu_*sqr(k_)/epsilon_;
 384
 385     #include "wallNonlinearViscosityI.H"
 386
 387     nonlinearStress_ = symm
 388     (
 389         pow(k_, 3.0)/sqr(epsilon_)
 390        *(
 391             Ctau1_/fEta_
 392            *(
 393                 (gradU_ & gradU_)
 394               + (gradU_ & gradU_)().T()
 395             )
 396           + Ctau2_/fEta_*(gradU_ & gradU_.T())
 397           + Ctau3_/fEta_*(gradU_.T() & gradU_)
 398         )
 399     );
 400 }
 401
 402
 403 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 404
 405 } // End namespace RASModels
 406 } // End namespace incompressible
 407 } // End namespace Foam
 408
 409 // ************************************************************************* //