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 "PDRkEpsilon.H"
00027 #include "../../dragModels/PDRDragModel/PDRDragModel.H"
00028 #include <OpenFOAM/addToRunTimeSelectionTable.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(PDRkEpsilon, 0);
00044 addToRunTimeSelectionTable(RASModel, PDRkEpsilon, dictionary);
00045
00046
00047
00048 PDRkEpsilon::PDRkEpsilon
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 C1_
00068 (
00069 dimensioned<scalar>::lookupOrAddToDict
00070 (
00071 "C1",
00072 coeffDict_,
00073 1.44
00074 )
00075 ),
00076 C2_
00077 (
00078 dimensioned<scalar>::lookupOrAddToDict
00079 (
00080 "C2",
00081 coeffDict_,
00082 1.92
00083 )
00084 ),
00085 sigmak_
00086 (
00087 dimensioned<scalar>::lookupOrAddToDict
00088 (
00089 "sigmak",
00090 coeffDict_,
00091 1.0
00092 )
00093 ),
00094 sigmaEps_
00095 (
00096 dimensioned<scalar>::lookupOrAddToDict
00097 (
00098 "sigmaEps",
00099 coeffDict_,
00100 1.3
00101 )
00102 ),
00103 Prt_
00104 (
00105 dimensioned<scalar>::lookupOrAddToDict
00106 (
00107 "Prt",
00108 coeffDict_,
00109 1.0
00110 )
00111 ),
00112
00113 k_
00114 (
00115 IOobject
00116 (
00117 "k",
00118 runTime_.timeName(),
00119 mesh_,
00120 IOobject::NO_READ,
00121 IOobject::AUTO_WRITE
00122 ),
00123 autoCreateK("k", mesh_)
00124 ),
00125
00126 epsilon_
00127 (
00128 IOobject
00129 (
00130 "epsilon",
00131 runTime_.timeName(),
00132 mesh_,
00133 IOobject::NO_READ,
00134 IOobject::AUTO_WRITE
00135 ),
00136 autoCreateEpsilon("epsilon", mesh_)
00137 ),
00138
00139 mut_
00140 (
00141 IOobject
00142 (
00143 "mut",
00144 runTime_.timeName(),
00145 mesh_,
00146 IOobject::NO_READ,
00147 IOobject::AUTO_WRITE
00148 ),
00149 autoCreateMut("mut", mesh_)
00150 ),
00151
00152 alphat_
00153 (
00154 IOobject
00155 (
00156 "alphat",
00157 runTime_.timeName(),
00158 mesh_,
00159 IOobject::NO_READ,
00160 IOobject::AUTO_WRITE
00161 ),
00162 autoCreateAlphat("alphat", mesh_)
00163 )
00164 {
00165 mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
00166 mut_.correctBoundaryConditions();
00167
00168 alphat_ = mut_/Prt_;
00169 alphat_.correctBoundaryConditions();
00170
00171 printCoeffs();
00172 }
00173
00174
00175
00176
00177 tmp<volSymmTensorField> PDRkEpsilon::R() const
00178 {
00179 return tmp<volSymmTensorField>
00180 (
00181 new volSymmTensorField
00182 (
00183 IOobject
00184 (
00185 "R",
00186 runTime_.timeName(),
00187 mesh_,
00188 IOobject::NO_READ,
00189 IOobject::NO_WRITE
00190 ),
00191 ((2.0/3.0)*I)*k_ - (mut_/rho_)*dev(twoSymm(fvc::grad(U_))),
00192 k_.boundaryField().types()
00193 )
00194 );
00195 }
00196
00197
00198 tmp<volSymmTensorField> PDRkEpsilon::devRhoReff() const
00199 {
00200 return tmp<volSymmTensorField>
00201 (
00202 new volSymmTensorField
00203 (
00204 IOobject
00205 (
00206 "devRhoReff",
00207 runTime_.timeName(),
00208 mesh_,
00209 IOobject::NO_READ,
00210 IOobject::NO_WRITE
00211 ),
00212 -muEff()*dev(twoSymm(fvc::grad(U_)))
00213 )
00214 );
00215 }
00216
00217
00218 tmp<fvVectorMatrix> PDRkEpsilon::divDevRhoReff(volVectorField& U) const
00219 {
00220 return
00221 (
00222 - fvm::laplacian(muEff(), U) - fvc::div(muEff()*dev2(fvc::grad(U)().T()))
00223 );
00224 }
00225
00226
00227 bool PDRkEpsilon::read()
00228 {
00229 if (RASModel::read())
00230 {
00231 Cmu_.readIfPresent(coeffDict_);
00232 C1_.readIfPresent(coeffDict_);
00233 C2_.readIfPresent(coeffDict_);
00234 sigmak_.readIfPresent(coeffDict());
00235 sigmaEps_.readIfPresent(coeffDict());
00236 Prt_.readIfPresent(coeffDict());
00237
00238 return true;
00239 }
00240 else
00241 {
00242 return false;
00243 }
00244 }
00245
00246
00247 void PDRkEpsilon::correct()
00248 {
00249 if (!turbulence_)
00250 {
00251
00252 mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
00253 mut_.correctBoundaryConditions();
00254
00255
00256 alphat_ = mut_/Prt_;
00257 alphat_.correctBoundaryConditions();
00258
00259 return;
00260 }
00261
00262 RASModel::correct();
00263
00264 volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));
00265
00266 if (mesh_.moving())
00267 {
00268 divU += fvc::div(mesh_.phi());
00269 }
00270
00271 tmp<volTensorField> tgradU = fvc::grad(U_);
00272 volScalarField G("RASModel::G", mut_*(tgradU() && dev(twoSymm(tgradU()))));
00273 tgradU.clear();
00274
00275
00276 epsilon_.boundaryField().updateCoeffs();
00277
00278
00279
00280 const volScalarField& betav = U_.db().lookupObject<volScalarField>("betav");
00281
00282 const PDRDragModel& drag =
00283 U_.db().lookupObject<PDRDragModel>("PDRDragModel");
00284
00285 volScalarField GR = drag.Gk();
00286
00287
00288 tmp<fvScalarMatrix> epsEqn
00289 (
00290 betav*fvm::ddt(rho_, epsilon_)
00291 + fvm::div(phi_, epsilon_)
00292 - fvm::laplacian(DepsilonEff(), epsilon_)
00293 ==
00294 C1_*(betav*G + GR)*epsilon_/k_
00295 - fvm::SuSp(((2.0/3.0)*C1_)*betav*rho_*divU, epsilon_)
00296 - fvm::Sp(C2_*betav*rho_*epsilon_/k_, epsilon_)
00297 );
00298
00299 epsEqn().relax();
00300
00301 epsEqn().boundaryManipulate(epsilon_.boundaryField());
00302
00303 solve(epsEqn);
00304 bound(epsilon_, epsilon0_);
00305
00306
00307
00308
00309 tmp<fvScalarMatrix> kEqn
00310 (
00311 betav*fvm::ddt(rho_, k_)
00312 + fvm::div(phi_, k_)
00313 - fvm::laplacian(DkEff(), k_)
00314 ==
00315 betav*G + GR
00316 - fvm::SuSp((2.0/3.0)*betav*rho_*divU, k_)
00317 - fvm::Sp(betav*rho_*epsilon_/k_, k_)
00318 );
00319
00320 kEqn().relax();
00321 solve(kEqn);
00322 bound(k_, k0_);
00323
00324
00325 mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
00326 mut_.correctBoundaryConditions();
00327
00328
00329 alphat_ = mut_/Prt_;
00330 alphat_.correctBoundaryConditions();
00331 }
00332
00333
00334
00335
00336 }
00337 }
00338 }
00339
00340