1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
28 #include <OpenFOAM/addToRunTimeSelectionTable.H>
29 #include <finiteVolume/wallFvPatch.H>
31 #include <compressibleRASModels/backwardsCompatibilityWallFunctions.H>
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 namespace compressible
42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44 defineTypeNameAndDebug(LRR, 0);
45 addToRunTimeSelectionTable(RASModel, LRR, dictionary);
47 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
51 const volScalarField& rho,
52 const volVectorField& U,
53 const surfaceScalarField& phi,
54 const basicThermo& thermophysicalModel
57 RASModel(typeName, rho, U, phi, thermophysicalModel),
61 dimensioned<scalar>::lookupOrAddToDict
70 dimensioned<scalar>::lookupOrAddToDict
79 dimensioned<scalar>::lookupOrAddToDict
88 dimensioned<scalar>::lookupOrAddToDict
97 dimensioned<scalar>::lookupOrAddToDict
106 dimensioned<scalar>::lookupOrAddToDict
115 dimensioned<scalar>::lookupOrAddToDict
124 dimensioned<scalar>::lookupOrAddToDict
133 dimensioned<scalar>::lookupOrAddToDict
142 dimensioned<scalar>::lookupOrAddToDict
151 dimensioned<scalar>::lookupOrAddToDict
169 autoCreateR("R", mesh_)
181 autoCreateK("k", mesh_)
193 autoCreateEpsilon("epsilon", mesh_)
205 autoCreateMut("mut", mesh_)
217 autoCreateAlphat("alphat", mesh_)
220 if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
225 "( const volScalarField&, const volVectorField&"
226 ", const surfaceScalarField&, incompressibleTransportModel&)"
227 ) << "couplingFactor = " << couplingFactor_
228 << " is not in range 0 - 1" << nl
232 mut_ = Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
233 mut_.correctBoundaryConditions();
236 alphat_.correctBoundaryConditions();
242 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
244 tmp<volSymmTensorField> LRR::devRhoReff() const
246 return tmp<volSymmTensorField>
248 new volSymmTensorField
258 rho_*R_ - mu()*dev(twoSymm(fvc::grad(U_)))
264 tmp<fvVectorMatrix> LRR::divDevRhoReff(volVectorField& U) const
266 if (couplingFactor_.value() > 0.0)
270 fvc::div(rho_*R_ + couplingFactor_*mut_*fvc::grad(U))
271 + fvc::laplacian((1.0 - couplingFactor_)*mut_, U)
272 - fvm::laplacian(muEff(), U)
273 - fvc::div(mu()*dev2(fvc::grad(U)().T()))
281 + fvc::laplacian(mut_, U)
282 - fvm::laplacian(muEff(), U)
283 - fvc::div(mu()*dev2(fvc::grad(U)().T()))
291 if (RASModel::read())
293 Cmu_.readIfPresent(coeffDict());
294 Clrr1_.readIfPresent(coeffDict());
295 Clrr2_.readIfPresent(coeffDict());
296 C1_.readIfPresent(coeffDict());
297 C2_.readIfPresent(coeffDict());
298 Cs_.readIfPresent(coeffDict());
299 Ceps_.readIfPresent(coeffDict());
300 sigmaR_.readIfPresent(coeffDict());
301 sigmaEps_.readIfPresent(coeffDict());
302 Prt_.readIfPresent(coeffDict());
303 couplingFactor_.readIfPresent(coeffDict());
305 if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
307 FatalErrorIn("LRR::read()")
308 << "couplingFactor = " << couplingFactor_
309 << " is not in range 0 - 1" << nl
326 // Re-calculate viscosity
327 mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
328 mut_.correctBoundaryConditions();
330 // Re-calculate thermal diffusivity
332 alphat_.correctBoundaryConditions();
339 volSymmTensorField P = -twoSymm(R_ & fvc::grad(U_));
340 volScalarField G("RASModel::G", 0.5*mag(tr(P)));
342 // Update espsilon and G at the wall
343 epsilon_.boundaryField().updateCoeffs();
345 // Dissipation equation
346 tmp<fvScalarMatrix> epsEqn
348 fvm::ddt(rho_, epsilon_)
349 + fvm::div(phi_, epsilon_)
350 //- fvm::laplacian(Ceps*rho_*(k_/epsilon_)*R_, epsilon_)
351 - fvm::laplacian(DepsilonEff(), epsilon_)
353 C1_*rho_*G*epsilon_/k_
354 - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)
359 epsEqn().boundaryManipulate(epsilon_.boundaryField());
362 bound(epsilon_, epsilon0_);
365 // Reynolds stress equation
367 const fvPatchList& patches = mesh_.boundary();
369 forAll(patches, patchi)
371 const fvPatch& curPatch = patches[patchi];
373 if (typeid(curPatch) == typeid(wallFvPatch))
375 forAll(curPatch, facei)
377 label faceCelli = curPatch.faceCells()[facei];
379 *= min(G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL), 100.0);
385 tmp<fvSymmTensorMatrix> REqn
389 //- fvm::laplacian(Cs*rho_*(k_/epsilon_)*R_, R_)
390 - fvm::laplacian(DREff(), R_)
391 + fvm::Sp(Clrr1_*rho_*epsilon_/k_, R_)
394 - (2.0/3.0*(1 - Clrr1_)*I)*rho_*epsilon_
403 dimensionedSymmTensor
409 k0_.value(), -GREAT, -GREAT,
420 // Re-calculate viscosity
421 mut_ = rho_*Cmu_*sqr(k_)/epsilon_;
422 mut_.correctBoundaryConditions();
424 // Re-calculate thermal diffusivity
426 alphat_.correctBoundaryConditions();
429 // Correct wall shear stresses
431 forAll(patches, patchi)
433 const fvPatch& curPatch = patches[patchi];
435 if (typeid(curPatch) == typeid(wallFvPatch))
437 symmTensorField& Rw = R_.boundaryField()[patchi];
439 const scalarField& rhow = rho_.boundaryField()[patchi];
440 const scalarField& mutw = mut_.boundaryField()[patchi];
442 vectorField snGradU = U_.boundaryField()[patchi].snGrad();
444 const vectorField& faceAreas
445 = mesh_.Sf().boundaryField()[patchi];
447 const scalarField& magFaceAreas
448 = mesh_.magSf().boundaryField()[patchi];
450 forAll(curPatch, facei)
452 // Calculate near-wall velocity gradient
454 = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
456 // Calculate near-wall shear-stress tensor
457 tensor tauw = -(mutw[facei]/rhow[facei])*2*dev(symm(gradUw));
459 // Reset the shear components of the stress tensor
460 Rw[facei].xy() = tauw.xy();
461 Rw[facei].xz() = tauw.xz();
462 Rw[facei].yz() = tauw.yz();
469 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
471 } // End namespace RASModels
472 } // End namespace compressible
473 } // End namespace Foam
475 // ************************ vim: set sw=4 sts=4 et: ************************ //