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 "SpalartAllmaras.H"
00027 #include <OpenFOAM/addToRunTimeSelectionTable.H>
00028
00029 #include <compressibleRASModels/backwardsCompatibilityWallFunctions.H>
00030
00031
00032
00033 namespace Foam
00034 {
00035 namespace compressible
00036 {
00037 namespace RASModels
00038 {
00039
00040
00041
00042 defineTypeNameAndDebug(SpalartAllmaras, 0);
00043 addToRunTimeSelectionTable(RASModel, SpalartAllmaras, dictionary);
00044
00045
00046
00047 tmp<volScalarField> SpalartAllmaras::chi() const
00048 {
00049 return rho_*nuTilda_/mu();
00050 }
00051
00052
00053 tmp<volScalarField> SpalartAllmaras::fv1(const volScalarField& chi) const
00054 {
00055 volScalarField chi3 = pow3(chi);
00056 return chi3/(chi3 + pow3(Cv1_));
00057 }
00058
00059
00060 tmp<volScalarField> SpalartAllmaras::fv2
00061 (
00062 const volScalarField& chi,
00063 const volScalarField& fv1
00064 ) const
00065 {
00066 return 1.0/pow3(scalar(1) + chi/Cv2_);
00067 }
00068
00069
00070 tmp<volScalarField> SpalartAllmaras::fv3
00071 (
00072 const volScalarField& chi,
00073 const volScalarField& fv1
00074 ) const
00075 {
00076 volScalarField chiByCv2 = (1/Cv2_)*chi;
00077
00078 return
00079 (scalar(1) + chi*fv1)
00080 *(1/Cv2_)
00081 *(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
00082 /pow3(scalar(1) + chiByCv2);
00083 }
00084
00085
00086 tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& Stilda) const
00087 {
00088 volScalarField r = min
00089 (
00090 nuTilda_
00091 /(
00092 max(Stilda, dimensionedScalar("SMALL", Stilda.dimensions(), SMALL))
00093 *sqr(kappa_*d_)
00094 ),
00095 scalar(10.0)
00096 );
00097 r.boundaryField() == 0.0;
00098
00099 volScalarField g = r + Cw2_*(pow6(r) - r);
00100
00101 return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
00102 }
00103
00104
00105
00106
00107 SpalartAllmaras::SpalartAllmaras
00108 (
00109 const volScalarField& rho,
00110 const volVectorField& U,
00111 const surfaceScalarField& phi,
00112 const basicThermo& thermophysicalModel
00113 )
00114 :
00115 RASModel(typeName, rho, U, phi, thermophysicalModel),
00116
00117 sigmaNut_
00118 (
00119 dimensioned<scalar>::lookupOrAddToDict
00120 (
00121 "sigmaNut",
00122 coeffDict_,
00123 0.66666
00124 )
00125 ),
00126 kappa_
00127 (
00128 dimensioned<scalar>::lookupOrAddToDict
00129 (
00130 "kappa",
00131 coeffDict_,
00132 0.41
00133 )
00134 ),
00135 Prt_
00136 (
00137 dimensioned<scalar>::lookupOrAddToDict
00138 (
00139 "Prt",
00140 coeffDict_,
00141 1.0
00142 )
00143 ),
00144
00145 Cb1_
00146 (
00147 dimensioned<scalar>::lookupOrAddToDict
00148 (
00149 "Cb1",
00150 coeffDict_,
00151 0.1355
00152 )
00153 ),
00154 Cb2_
00155 (
00156 dimensioned<scalar>::lookupOrAddToDict
00157 (
00158 "Cb2",
00159 coeffDict_,
00160 0.622
00161 )
00162 ),
00163 Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
00164 Cw2_
00165 (
00166 dimensioned<scalar>::lookupOrAddToDict
00167 (
00168 "Cw2",
00169 coeffDict_,
00170 0.3
00171 )
00172 ),
00173 Cw3_
00174 (
00175 dimensioned<scalar>::lookupOrAddToDict
00176 (
00177 "Cw3",
00178 coeffDict_,
00179 2.0
00180 )
00181 ),
00182 Cv1_
00183 (
00184 dimensioned<scalar>::lookupOrAddToDict
00185 (
00186 "Cv1",
00187 coeffDict_,
00188 7.1
00189 )
00190 ),
00191 Cv2_
00192 (
00193 dimensioned<scalar>::lookupOrAddToDict
00194 (
00195 "Cv2",
00196 coeffDict_,
00197 5.0
00198 )
00199 ),
00200
00201 nuTilda_
00202 (
00203 IOobject
00204 (
00205 "nuTilda",
00206 runTime_.timeName(),
00207 mesh_,
00208 IOobject::MUST_READ,
00209 IOobject::AUTO_WRITE
00210 ),
00211 mesh_
00212 ),
00213
00214 mut_
00215 (
00216 IOobject
00217 (
00218 "mut",
00219 runTime_.timeName(),
00220 mesh_,
00221 IOobject::MUST_READ,
00222 IOobject::AUTO_WRITE
00223 ),
00224 mesh_
00225 ),
00226
00227 alphat_
00228 (
00229 IOobject
00230 (
00231 "alphat",
00232 runTime_.timeName(),
00233 mesh_,
00234 IOobject::NO_READ,
00235 IOobject::AUTO_WRITE
00236 ),
00237 autoCreateAlphat("alphat", mesh_)
00238 ),
00239
00240 d_(mesh_)
00241 {
00242 alphat_ = mut_/Prt_;
00243 alphat_.correctBoundaryConditions();
00244
00245 printCoeffs();
00246 }
00247
00248
00249
00250
00251 tmp<volSymmTensorField> SpalartAllmaras::R() const
00252 {
00253 return tmp<volSymmTensorField>
00254 (
00255 new volSymmTensorField
00256 (
00257 IOobject
00258 (
00259 "R",
00260 runTime_.timeName(),
00261 mesh_,
00262 IOobject::NO_READ,
00263 IOobject::NO_WRITE
00264 ),
00265 ((2.0/3.0)*I)*k() - (mut_/rho_)*dev(twoSymm(fvc::grad(U_)))
00266 )
00267 );
00268 }
00269
00270
00271 tmp<volSymmTensorField> SpalartAllmaras::devRhoReff() const
00272 {
00273 return tmp<volSymmTensorField>
00274 (
00275 new volSymmTensorField
00276 (
00277 IOobject
00278 (
00279 "devRhoReff",
00280 runTime_.timeName(),
00281 mesh_,
00282 IOobject::NO_READ,
00283 IOobject::NO_WRITE
00284 ),
00285 -muEff()*dev(twoSymm(fvc::grad(U_)))
00286 )
00287 );
00288 }
00289
00290
00291 tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff(volVectorField& U) const
00292 {
00293 volScalarField muEff_ = muEff();
00294
00295 return
00296 (
00297 - fvm::laplacian(muEff_, U)
00298 - fvc::div(muEff_*dev2(fvc::grad(U)().T()))
00299 );
00300 }
00301
00302
00303 bool SpalartAllmaras::read()
00304 {
00305 if (RASModel::read())
00306 {
00307 sigmaNut_.readIfPresent(coeffDict());
00308 kappa_.readIfPresent(coeffDict());
00309 Prt_.readIfPresent(coeffDict());
00310
00311 Cb1_.readIfPresent(coeffDict());
00312 Cb2_.readIfPresent(coeffDict());
00313 Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
00314 Cw2_.readIfPresent(coeffDict());
00315 Cw3_.readIfPresent(coeffDict());
00316 Cv1_.readIfPresent(coeffDict());
00317 Cv2_.readIfPresent(coeffDict());
00318
00319 return true;
00320 }
00321 else
00322 {
00323 return false;
00324 }
00325 }
00326
00327
00328 void SpalartAllmaras::correct()
00329 {
00330 if (!turbulence_)
00331 {
00332
00333 mut_ = rho_*nuTilda_*fv1(chi());
00334 mut_.correctBoundaryConditions();
00335
00336
00337 alphat_ = mut_/Prt_;
00338 alphat_.correctBoundaryConditions();
00339
00340 return;
00341 }
00342
00343 RASModel::correct();
00344
00345 if (mesh_.changing())
00346 {
00347 d_.correct();
00348 }
00349
00350 volScalarField chi = this->chi();
00351 volScalarField fv1 = this->fv1(chi);
00352
00353 volScalarField Stilda =
00354 fv3(chi, fv1)*::sqrt(2.0)*mag(skew(fvc::grad(U_)))
00355 + fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_);
00356
00357 tmp<fvScalarMatrix> nuTildaEqn
00358 (
00359 fvm::ddt(rho_, nuTilda_)
00360 + fvm::div(phi_, nuTilda_)
00361 - fvm::laplacian(DnuTildaEff(), nuTilda_)
00362 - Cb2_/sigmaNut_*rho_*magSqr(fvc::grad(nuTilda_))
00363 ==
00364 Cb1_*rho_*Stilda*nuTilda_
00365 - fvm::Sp(Cw1_*fw(Stilda)*nuTilda_*rho_/sqr(d_), nuTilda_)
00366 );
00367
00368 nuTildaEqn().relax();
00369 solve(nuTildaEqn);
00370 bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
00371 nuTilda_.correctBoundaryConditions();
00372
00373
00374 mut_.internalField() = fv1*nuTilda_.internalField()*rho_.internalField();
00375 mut_.correctBoundaryConditions();
00376
00377
00378 alphat_ = mut_/Prt_;
00379 alphat_.correctBoundaryConditions();
00380 }
00381
00382
00383
00384
00385 }
00386 }
00387 }
00388
00389