00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
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
00042
00043 defineTypeNameAndDebug(LRR, 0);
00044 addToRunTimeSelectionTable(RASModel, LRR, dictionary);
00045
00046
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
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
00326 mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
00327 mut_.correctBoundaryConditions();
00328
00329
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
00342 epsilon_.boundaryField().updateCoeffs();
00343
00344
00345 tmp<fvScalarMatrix> epsEqn
00346 (
00347 fvm::ddt(rho_, epsilon_)
00348 + fvm::div(phi_, epsilon_)
00349
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
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
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
00420 mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
00421 mut_.correctBoundaryConditions();
00422
00423
00424 alphat_ = mut_/Prt_;
00425 alphat_.correctBoundaryConditions();
00426
00427
00428
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
00452 tensor gradUw
00453 = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
00454
00455
00456 tensor tauw = -(mutw[facei]/rhow[facei])*2*dev(symm(gradUw));
00457
00458
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 }
00471 }
00472 }
00473
00474