FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

LRR.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "LRR.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028 #include <finiteVolume/wallFvPatch.H>
00029 
00030 #include <compressibleRASModels/backwardsCompatibilityWallFunctions.H>
00031 
00032 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00033 
00034 namespace Foam
00035 {
00036 namespace compressible
00037 {
00038 namespace RASModels
00039 {
00040 
00041 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
00042 
00043 defineTypeNameAndDebug(LRR, 0);
00044 addToRunTimeSelectionTable(RASModel, LRR, dictionary);
00045 
00046 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00047 
00048 LRR::LRR
00049 (
00050     const volScalarField& rho,
00051     const volVectorField& U,
00052     const surfaceScalarField& phi,
00053     const basicThermo& thermophysicalModel
00054 )
00055 :
00056     RASModel(typeName, rho, U, phi, thermophysicalModel),
00057 
00058     Cmu_
00059     (
00060         dimensioned<scalar>::lookupOrAddToDict
00061         (
00062             "Cmu",
00063             coeffDict_,
00064             0.09
00065         )
00066     ),
00067     Clrr1_
00068     (
00069         dimensioned<scalar>::lookupOrAddToDict
00070         (
00071             "Clrr1",
00072             coeffDict_,
00073             1.8
00074         )
00075     ),
00076     Clrr2_
00077     (
00078         dimensioned<scalar>::lookupOrAddToDict
00079         (
00080             "Clrr2",
00081             coeffDict_,
00082             0.6
00083         )
00084     ),
00085     C1_
00086     (
00087         dimensioned<scalar>::lookupOrAddToDict
00088         (
00089             "C1",
00090             coeffDict_,
00091             1.44
00092         )
00093     ),
00094     C2_
00095     (
00096         dimensioned<scalar>::lookupOrAddToDict
00097         (
00098             "C2",
00099             coeffDict_,
00100             1.92
00101         )
00102     ),
00103     Cs_
00104     (
00105         dimensioned<scalar>::lookupOrAddToDict
00106         (
00107             "Cs",
00108             coeffDict_,
00109             0.25
00110         )
00111     ),
00112     Ceps_
00113     (
00114         dimensioned<scalar>::lookupOrAddToDict
00115         (
00116             "Ceps",
00117             coeffDict_,
00118             0.15
00119         )
00120     ),
00121     couplingFactor_
00122     (
00123         dimensioned<scalar>::lookupOrAddToDict
00124         (
00125             "couplingFactor",
00126             coeffDict_,
00127             0.0
00128         )
00129     ),
00130     sigmaR_
00131     (
00132         dimensioned<scalar>::lookupOrAddToDict
00133         (
00134             "sigmaR",
00135             coeffDict_,
00136             0.81967
00137         )
00138     ),
00139     sigmaEps_
00140     (
00141         dimensioned<scalar>::lookupOrAddToDict
00142         (
00143             "sigmaEps",
00144             coeffDict_,
00145             1.3
00146         )
00147     ),
00148     Prt_
00149     (
00150         dimensioned<scalar>::lookupOrAddToDict
00151         (
00152             "Prt",
00153             coeffDict_,
00154             1.0
00155         )
00156     ),
00157 
00158     R_
00159     (
00160         IOobject
00161         (
00162             "R",
00163             runTime_.timeName(),
00164             mesh_,
00165             IOobject::MUST_READ,
00166             IOobject::AUTO_WRITE
00167         ),
00168         autoCreateR("R", mesh_)
00169     ),
00170     k_
00171     (
00172         IOobject
00173         (
00174             "k",
00175             runTime_.timeName(),
00176             mesh_,
00177             IOobject::NO_READ,
00178             IOobject::AUTO_WRITE
00179         ),
00180         autoCreateK("k", mesh_)
00181     ),
00182     epsilon_
00183     (
00184         IOobject
00185         (
00186             "epsilon",
00187             runTime_.timeName(),
00188             mesh_,
00189             IOobject::NO_READ,
00190             IOobject::AUTO_WRITE
00191         ),
00192         autoCreateEpsilon("epsilon", mesh_)
00193     ),
00194     mut_
00195     (
00196         IOobject
00197         (
00198             "mut",
00199             runTime_.timeName(),
00200             mesh_,
00201             IOobject::NO_READ,
00202             IOobject::AUTO_WRITE
00203         ),
00204         autoCreateMut("mut", mesh_)
00205     ),
00206     alphat_
00207     (
00208         IOobject
00209         (
00210             "alphat",
00211             runTime_.timeName(),
00212             mesh_,
00213             IOobject::NO_READ,
00214             IOobject::AUTO_WRITE
00215         ),
00216         autoCreateAlphat("alphat", mesh_)
00217     )
00218 {
00219     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
00220     {
00221         FatalErrorIn
00222         (
00223             "LRR::LRR"
00224             "( const volScalarField&, const volVectorField&"
00225             ", const surfaceScalarField&, incompressibleTransportModel&)"
00226         )   << "couplingFactor = " << couplingFactor_
00227             << " is not in range 0 - 1" << nl
00228             << exit(FatalError);
00229     }
00230 
00231     mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
00232     mut_.correctBoundaryConditions();
00233 
00234     alphat_ = mut_/Prt_;
00235     alphat_.correctBoundaryConditions();
00236 
00237     printCoeffs();
00238 }
00239 
00240 
00241 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00242 
00243 tmp<volSymmTensorField> LRR::devRhoReff() const
00244 {
00245     return tmp<volSymmTensorField>
00246     (
00247         new volSymmTensorField
00248         (
00249             IOobject
00250             (
00251                 "devRhoReff",
00252                 runTime_.timeName(),
00253                 mesh_,
00254                 IOobject::NO_READ,
00255                 IOobject::NO_WRITE
00256             ),
00257             rho_*R_ - mu()*dev(twoSymm(fvc::grad(U_)))
00258         )
00259     );
00260 }
00261 
00262 
00263 tmp<fvVectorMatrix> LRR::divDevRhoReff(volVectorField& U) const
00264 {
00265     if (couplingFactor_.value() > 0.0)
00266     {
00267         return
00268         (
00269             fvc::div(rho_*R_ + couplingFactor_*mut_*fvc::grad(U))
00270           + fvc::laplacian((1.0 - couplingFactor_)*mut_, U)
00271           - fvm::laplacian(muEff(), U)
00272           - fvc::div(mu()*dev2(fvc::grad(U)().T()))
00273         );
00274     }
00275     else
00276     {
00277         return
00278         (
00279             fvc::div(rho_*R_)
00280           + fvc::laplacian(mut_, U)
00281           - fvm::laplacian(muEff(), U)
00282           - fvc::div(mu()*dev2(fvc::grad(U)().T()))
00283         );
00284     }
00285 }
00286 
00287 
00288 bool LRR::read()
00289 {
00290     if (RASModel::read())
00291     {
00292         Cmu_.readIfPresent(coeffDict());
00293         Clrr1_.readIfPresent(coeffDict());
00294         Clrr2_.readIfPresent(coeffDict());
00295         C1_.readIfPresent(coeffDict());
00296         C2_.readIfPresent(coeffDict());
00297         Cs_.readIfPresent(coeffDict());
00298         Ceps_.readIfPresent(coeffDict());
00299         sigmaR_.readIfPresent(coeffDict());
00300         sigmaEps_.readIfPresent(coeffDict());
00301         Prt_.readIfPresent(coeffDict());
00302         couplingFactor_.readIfPresent(coeffDict());
00303 
00304         if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
00305         {
00306             FatalErrorIn("LRR::read()")
00307                 << "couplingFactor = " << couplingFactor_
00308                 << " is not in range 0 - 1" << nl
00309                 << exit(FatalError);
00310         }
00311 
00312         return true;
00313     }
00314     else
00315     {
00316         return false;
00317     }
00318 }
00319 
00320 
00321 void LRR::correct()
00322 {
00323     if (!turbulence_)
00324     {
00325         // Re-calculate viscosity
00326         mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
00327         mut_.correctBoundaryConditions();
00328 
00329         // Re-calculate thermal diffusivity
00330         alphat_ = mut_/Prt_;
00331         alphat_.correctBoundaryConditions();
00332 
00333         return;
00334     }
00335 
00336     RASModel::correct();
00337 
00338     volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
00339     volScalarField G("RASModel::G", 0.5*mag(tr(P)));
00340 
00341     // Update espsilon and G at the wall
00342     epsilon_.boundaryField().updateCoeffs();
00343 
00344     // Dissipation equation
00345     tmp<fvScalarMatrix> epsEqn
00346     (
00347         fvm::ddt(rho_, epsilon_)
00348       + fvm::div(phi_, epsilon_)
00349     //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
00350       - fvm::laplacian(DepsilonEff(), epsilon_)
00351      ==
00352         C1_*rho_*G*epsilon_/k_
00353       - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
00354     );
00355 
00356     epsEqn().relax();
00357 
00358     epsEqn().boundaryManipulate(epsilon_.boundaryField());
00359 
00360     solve(epsEqn);
00361     bound(epsilon_, epsilon0_);
00362 
00363 
00364     // Reynolds stress equation
00365 
00366     const fvPatchList& patches = mesh_.boundary();
00367 
00368     forAll(patches, patchi)
00369     {
00370         const fvPatch& curPatch = patches[patchi];
00371 
00372         if (isA<wallFvPatch>(curPatch))
00373         {
00374             forAll(curPatch, facei)
00375             {
00376                 label faceCelli = curPatch.faceCells()[facei];
00377                 P[faceCelli]
00378                     *= min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 100.0);
00379             }
00380         }
00381     }
00382 
00383 
00384     tmp<fvSymmTensorMatrix> REqn
00385     (
00386         fvm::ddt(rho_, R_)
00387       + fvm::div(phi_, R_)
00388     //- fvm::laplacian(Cs*rho_*(k_/epsilon_)*R_, R_)
00389       - fvm::laplacian(DREff(), R_)
00390       + fvm::Sp(Clrr1_*rho_*epsilon_/k_, R_)
00391      ==
00392         rho_*P
00393       - (2.0/3.0*(1 - Clrr1_)*I)*rho_*epsilon_
00394       - Clrr2_*rho_*dev(P)
00395     );
00396 
00397     REqn().relax();
00398     solve(REqn);
00399 
00400     R_.max
00401     (
00402         dimensionedSymmTensor
00403         (
00404             "zero",
00405             R_.dimensions(),
00406             symmTensor
00407             (
00408                 k0_.value(), -GREAT, -GREAT,
00409                 k0_.value(), -GREAT,
00410                 k0_.value()
00411             )
00412         )
00413     );
00414 
00415     k_ = 0.5*tr(R_);
00416     bound(k_, k0_);
00417 
00418 
00419     // Re-calculate viscosity
00420     mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
00421     mut_.correctBoundaryConditions();
00422 
00423     // Re-calculate thermal diffusivity
00424     alphat_ = mut_/Prt_;
00425     alphat_.correctBoundaryConditions();
00426 
00427 
00428     // Correct wall shear stresses
00429 
00430     forAll(patches, patchi)
00431     {
00432         const fvPatch& curPatch = patches[patchi];
00433 
00434         if (isA<wallFvPatch>(curPatch))
00435         {
00436             symmTensorField& Rw = R_.boundaryField()[patchi];
00437 
00438             const scalarField& rhow = rho_.boundaryField()[patchi];
00439             const scalarField& mutw = mut_.boundaryField()[patchi];
00440 
00441             vectorField snGradU = U_.boundaryField()[patchi].snGrad();
00442 
00443             const vectorField& faceAreas
00444                 = mesh_.Sf().boundaryField()[patchi];
00445 
00446             const scalarField& magFaceAreas
00447                 = mesh_.magSf().boundaryField()[patchi];
00448 
00449             forAll(curPatch, facei)
00450             {
00451                 // Calculate near-wall velocity gradient
00452                 tensor gradUw
00453                     = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
00454 
00455                 // Calculate near-wall shear-stress tensor
00456                 tensor tauw = -(mutw[facei]/rhow[facei])*2*dev(symm(gradUw));
00457 
00458                 // Reset the shear components of the stress tensor
00459                 Rw[facei].xy() = tauw.xy();
00460                 Rw[facei].xz() = tauw.xz();
00461                 Rw[facei].yz() = tauw.yz();
00462             }
00463         }
00464     }
00465 }
00466 
00467 
00468 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00469 
00470 } // End namespace RASModels
00471 } // End namespace compressible
00472 } // End namespace Foam
00473 
00474 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines