zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Removed turbulence model coefficients from the dictionaries to allow them to default.

Browse Source
This commit is contained in:
henry
2009-07-21 19:05:30 +01:00
parent fb9e5bc4fe
commit b3d455de80
117 changed files with 331 additions and 6167 deletions
Download Patch File Download Diff File Expand all files Collapse all files

79
applications/solvers/multiphase/bubbleFoam/createRASTurbulence.H
View File

@ -16,52 +16,93 @@
RASProperties.lookup("turbulence")
);
dictionary kEpsilonCoeffs
dictionary kEpsilonDict
(
RASProperties.subDict("kEpsilonCoeffs")
RASProperties.subDictPtr("kEpsilonCoeffs")
);
scalar Cmu
dimensionedScalar Cmu
(
readScalar(kEpsilonCoeffs.lookup("Cmu"))
dimensionedScalar::lookupOrAddToDict
(
"Cmu",
kEpsilonDict,
0.09
)
);
scalar C1
dimensionedScalar C1
(
readScalar(kEpsilonCoeffs.lookup("C1"))
dimensionedScalar::lookupOrAddToDict
(
"C1",
kEpsilonDict,
1.44
)
);
scalar C2
dimensionedScalar C2
(
readScalar(kEpsilonCoeffs.lookup("C2"))
dimensionedScalar::lookupOrAddToDict
(
"C2",
kEpsilonDict,
1.92
)
);
scalar alphak
dimensionedScalar alphak
(
readScalar(kEpsilonCoeffs.lookup("alphak"))
dimensionedScalar::lookupOrAddToDict
(
"alphaEps",
kEpsilonDict,
1.0
)
);
scalar alphaEps
dimensionedScalar alphaEps
(
readScalar(kEpsilonCoeffs.lookup("alphaEps"))
dimensionedScalar::lookupOrAddToDict
(
"alphaEps",
kEpsilonDict,
0.76923
)
);
dictionary wallFunctionCoeffs
dictionary wallFunctionDict
(
RASProperties.subDict("wallFunctionCoeffs")
RASProperties.subDictPtr("wallFunctionCoeffs")
);
scalar kappa
dimensionedScalar kappa
(
readScalar(wallFunctionCoeffs.lookup("kappa"))
dimensionedScalar::lookupOrAddToDict
(
"kappa",
wallFunctionDict,
0.4187
)
);
scalar E
dimensionedScalar E
(
readScalar(wallFunctionCoeffs.lookup("E"))
dimensionedScalar::lookupOrAddToDict
(
"E",
wallFunctionDict,
9.0
)
);
if (RASProperties.lookupOrDefault<Switch>("printCoeffs", false))
{
Info<< "kEpsilonCoeffs" << kEpsilonDict << nl
<< "wallFunctionCoeffs" << wallFunctionDict << endl;
}
nearWallDist y(mesh);

12
applications/solvers/multiphase/bubbleFoam/kEpsilon.H
View File

@ -16,7 +16,11 @@ if(turbulence)
(
fvm::ddt(beta, epsilon)
+ fvm::div(phib, epsilon)
- fvm::laplacian(alphaEps*nuEffb, epsilon)
- fvm::laplacian
(
alphaEps*nuEffb, epsilon,
"laplacian(DepsilonEff,epsilon)"
)
==
C1*beta*G*epsilon/k
- fvm::Sp(C2*beta*epsilon/k, epsilon)
@ -35,7 +39,11 @@ if(turbulence)
(
fvm::ddt(beta, k)
+ fvm::div(phib, k)
- fvm::laplacian(alphak*nuEffb, k)
- fvm::laplacian
(
alphak*nuEffb, k,
"laplacian(DkEff,k)"
)
==
beta*G
- fvm::Sp(beta*epsilon/k, k)

9
applications/solvers/multiphase/bubbleFoam/wallFunctions.H
View File

@ -1,8 +1,9 @@
{
labelList cellBoundaryFaceCount(epsilon.size(), 0);
scalar Cmu25 = ::pow(Cmu, 0.25);
scalar Cmu75 = ::pow(Cmu, 0.75);
scalar Cmu25 = ::pow(Cmu.value(), 0.25);
scalar Cmu75 = ::pow(Cmu.value(), 0.75);
scalar kappa_ = kappa.value();
const fvPatchList& patches = mesh.boundary();
@ -53,14 +54,14 @@
epsilon[faceCelli] +=
Cmu75*::pow(k[faceCelli], 1.5)
/(kappa*y[patchi][facei]);
/(kappa_*y[patchi][facei]);
if (yPlus > 11.6)
{
G[faceCelli] +=
nuw[facei]*magFaceGradU[facei]
*Cmu25*::sqrt(k[faceCelli])
/(kappa*y[patchi][facei]);
/(kappa_*y[patchi][facei]);
}
}
}

13
applications/solvers/multiphase/bubbleFoam/wallViscosity.H
View File

@ -1,5 +1,8 @@
{
scalar Cmu25 = ::pow(Cmu, 0.25);
scalar Cmu25 = ::pow(Cmu.value(), 0.25);
scalar kappa_ = kappa.value();
scalar E_ = E.value();
scalar nub_ = nub.value();
const fvPatchList& patches = mesh.boundary();
@ -19,14 +22,14 @@
scalar yPlus =
Cmu25*y[patchi][facei]
*::sqrt(k[faceCelli])
/nub.value();
/nub_;
if (yPlus > 11.6)
{
nutw[facei] =
yPlus*nub.value()*kappa
/::log(E*yPlus)
- nub.value();
yPlus*nub_*kappa_
/::log(E_*yPlus)
- nub_;
}
else
{

88
applications/solvers/multiphase/settlingFoam/createFields.H
View File

@ -228,57 +228,103 @@
RASProperties.lookup("turbulence")
);
const dictionary& kEpsilonCoeffs
dictionary kEpsilonDict
(
RASProperties.subDict("kEpsilonCoeffs")
RASProperties.subDictPtr("kEpsilonCoeffs")
);
scalar Cmu
dimensionedScalar Cmu
(
readScalar(kEpsilonCoeffs.lookup("Cmu"))
dimensionedScalar::lookupOrAddToDict
(
"Cmu",
kEpsilonDict,
0.09
)
);
scalar C1
dimensionedScalar C1
(
readScalar(kEpsilonCoeffs.lookup("C1"))
dimensionedScalar::lookupOrAddToDict
(
"C1",
kEpsilonDict,
1.44
)
);
scalar C2
dimensionedScalar C2
(
readScalar(kEpsilonCoeffs.lookup("C2"))
dimensionedScalar::lookupOrAddToDict
(
"C2",
kEpsilonDict,
1.92
)
);
scalar C3
dimensionedScalar C3
(
readScalar(kEpsilonCoeffs.lookup("C3"))
dimensionedScalar::lookupOrAddToDict
(
"C3",
kEpsilonDict,
0.85
)
);
scalar alphak
dimensionedScalar alphak
(
readScalar(kEpsilonCoeffs.lookup("alphak"))
dimensionedScalar::lookupOrAddToDict
(
"alphaEps",
kEpsilonDict,
1.0
)
);
scalar alphaEps
dimensionedScalar alphaEps
(
readScalar(kEpsilonCoeffs.lookup("alphaEps"))
dimensionedScalar::lookupOrAddToDict
(
"alphaEps",
kEpsilonDict,
0.76923
)
);
const dictionary& wallFunctionCoeffs
dictionary wallFunctionDict
(
RASProperties.subDict("wallFunctionCoeffs")
RASProperties.subDictPtr("wallFunctionCoeffs")
);
scalar kappa
dimensionedScalar kappa
(
readScalar(wallFunctionCoeffs.lookup("kappa"))
dimensionedScalar::lookupOrAddToDict
(
"kappa",
wallFunctionDict,
0.4187
)
);
scalar E
dimensionedScalar E
(
readScalar(wallFunctionCoeffs.lookup("E"))
dimensionedScalar::lookupOrAddToDict
(
"E",
wallFunctionDict,
9.0
)
);
if (RASProperties.lookupOrDefault<Switch>("printCoeffs", false))
{
Info<< "kEpsilonCoeffs" << kEpsilonDict << nl
<< "wallFunctionCoeffs" << wallFunctionDict << endl;
}
nearWallDist y(mesh);
Info<< "Reading field k\n" << endl;

8
applications/solvers/multiphase/settlingFoam/kEpsilon.H
View File

@ -26,7 +26,7 @@ if(turbulence)
+ fvm::div(phi, epsilon)
- fvm::laplacian
(
alphaEps*mut + mul, epsilon,
alphaEps*mut + mul, epsilon,
"laplacian(DepsilonEff,epsilon)"
)
==
@ -49,7 +49,11 @@ if(turbulence)
(
fvm::ddt(rho, k)
+ fvm::div(phi, k)
- fvm::laplacian(alphak*mut + mul, k, "laplacian(DkEff,k)")
- fvm::laplacian
(
alphak*mut + mul, k,
"laplacian(DkEff,k)"
)
==
G
- fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k)

9
applications/solvers/multiphase/settlingFoam/wallFunctions.H
View File

@ -1,8 +1,9 @@
{
labelList cellBoundaryFaceCount(epsilon.size(), 0);
scalar Cmu25 = ::pow(Cmu, 0.25);
scalar Cmu75 = ::pow(Cmu, 0.75);
scalar Cmu25 = ::pow(Cmu.value(), 0.25);
scalar Cmu75 = ::pow(Cmu.value(), 0.75);
scalar kappa_ = kappa.value();
const fvPatchList& patches = mesh.boundary();
@ -55,14 +56,14 @@
epsilon[faceCelli] +=
Cmu75*rho[faceCelli]*::pow(k[faceCelli], 1.5)
/(kappa*y[patchi][facei]);
/(kappa_*y[patchi][facei]);
if (yPlus > 11.6)
{
G[faceCelli] +=
mutw[facei]*magFaceGradU[facei]
*Cmu25*::sqrt(k[faceCelli])
/(kappa*y[patchi][facei]);
/(kappa_*y[patchi][facei]);
}
}
}

6
applications/solvers/multiphase/settlingFoam/wallViscosity.H
View File

@ -1,5 +1,7 @@
{
scalar Cmu25 = ::pow(Cmu, 0.25);
scalar Cmu25 = ::pow(Cmu.value(), 0.25);
scalar kappa_ = kappa.value();
scalar E_ = E.value();
const fvPatchList& patches = mesh.boundary();
@ -26,7 +28,7 @@
{
mutw[facei] =
muw[facei]
*(yPlus*kappa/::log(E*yPlus) - 1);
*(yPlus*kappa_/::log(E_*yPlus) - 1);
}
else
{

3
applications/solvers/multiphase/twoPhaseEulerFoam/Make/options
View File

@ -1,4 +1,5 @@
EXE_INC = \
-I../bubbleFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-IturbulenceModel \
@ -6,7 +7,7 @@ EXE_INC = \
-IinterfacialModels/lnInclude \
-IphaseModel/lnInclude \
-Iaveraging
EXE_LIBS = \
-lEulerianInterfacialModels \
-lfiniteVolume \

137
applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
View File

@ -91,34 +91,6 @@
alpha*Ua + beta*Ub
);
Info<< "Reading field k\n" << endl;
volScalarField k
(
IOobject
(
"k",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field epsilon\n" << endl;
volScalarField epsilon
(
IOobject
(
"epsilon",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
dimensionedScalar Cvm
(
transportProperties.lookup("Cvm")
@ -156,114 +128,7 @@
alpha*rhoa + beta*rhob
);
IOdictionary RASProperties
(
IOobject
(
"RASProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
Switch turbulence
(
RASProperties.lookup("turbulence")
);
dictionary kEpsilonCoeffs
(
RASProperties.subDict("kEpsilonCoeffs")
);
scalar Cmu
(
readScalar(kEpsilonCoeffs.lookup("Cmu"))
);
scalar C1
(
readScalar(kEpsilonCoeffs.lookup("C1"))
);
scalar C2
(
readScalar(kEpsilonCoeffs.lookup("C2"))
);
scalar alphak
(
readScalar(kEpsilonCoeffs.lookup("alphak"))
);
scalar alphaEps
(
readScalar(kEpsilonCoeffs.lookup("alphaEps"))
);
dictionary wallFunctionCoeffs
(
RASProperties.subDict("wallFunctionCoeffs")
);
scalar kappa
(
readScalar(wallFunctionCoeffs.lookup("kappa"))
);
scalar E
(
readScalar(wallFunctionCoeffs.lookup("E"))
);
nearWallDist y(mesh);
Info<< "Calculating field nutb\n" << endl;
volScalarField nutb
(
IOobject
(
"nutb",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Cmu*sqr(k)/epsilon
);
Info<< "Calculating field nuEffa\n" << endl;
volScalarField nuEffa
(
IOobject
(
"nuEffa",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
sqr(Ct)*nutb + nua
);
Info<< "Calculating field nuEffb\n" << endl;
volScalarField nuEffb
(
IOobject
(
"nuEffb",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
nutb + nub
);
#include "createRASTurbulence.H"
Info<< "Calculating field DDtUa and DDtUb\n" << endl;

61
applications/solvers/multiphase/twoPhaseEulerFoam/turbulenceModel/kEpsilon.H
View File

@ -1,61 +0,0 @@
if(turbulence)
{
if (mesh.changing())
{
y.correct();
}
tmp<volTensorField> tgradUb = fvc::grad(Ub);
volScalarField G = 2*nutb*(tgradUb() && dev(symm(tgradUb())));
tgradUb.clear();
# include "wallFunctions.H"
// Dissipation equation
fvScalarMatrix epsEqn
(
fvm::ddt(beta, epsilon)
+ fvm::div(phib, epsilon)
- fvm::laplacian
(
alphaEps*nuEffb,
epsilon,
"laplacian((alphaEps*nuEffb),epsilon)"
)
==
C1*beta*G*epsilon/k
- fvm::Sp(C2*beta*epsilon/k, epsilon)
);
# include "wallDissipation.H"
epsEqn.relax();
epsEqn.solve();
epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15));
// Turbulent kinetic energy equation
fvScalarMatrix kEqn
(
fvm::ddt(beta, k)
+ fvm::div(phib, k)
- fvm::laplacian(alphak*nuEffb, k, "laplacian((alphak*nuEffb),k)")
==
beta*G
- fvm::Sp(beta*epsilon/k, k)
);
kEqn.relax();
kEqn.solve();
k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8));
//- Re-calculate turbulence viscosity
nutb = Cmu*sqr(k)/epsilon;
# include "wallViscosity.H"
}
nuEffa = sqr(Ct)*nutb + nua;
nuEffb = nutb + nub;

51
applications/solvers/multiphase/twoPhaseEulerFoam/turbulenceModel/wallDissipation.H
View File

@ -1,51 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Global
wallDissipation
Description
Set wall dissipation in the epsilon matrix
\*---------------------------------------------------------------------------*/
{
const fvPatchList& patches = mesh.boundary();
forAll(patches, patchi)
{
const fvPatch& p = patches[patchi];
if (isType<wallFvPatch>(p))
{
epsEqn.setValues
(
p.faceCells(),
epsilon.boundaryField()[patchi].patchInternalField()
);
}
}
}
// ************************************************************************* //

86
applications/solvers/multiphase/twoPhaseEulerFoam/turbulenceModel/wallFunctions.H
View File

@ -1,86 +0,0 @@
{
labelList cellBoundaryFaceCount(epsilon.size(), 0);
scalar Cmu25 = ::pow(Cmu, 0.25);
scalar Cmu75 = ::pow(Cmu, 0.75);
const fvPatchList& patches = mesh.boundary();
//- Initialise the near-wall P field to zero
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (isType<wallFvPatch>(currPatch))
{
forAll(currPatch, facei)
{
label faceCelli = currPatch.faceCells()[facei];
epsilon[faceCelli] = 0.0;
G[faceCelli] = 0.0;
}
}
}
//- Accumulate the wall face contributions to epsilon and G
// Increment cellBoundaryFaceCount for each face for averaging
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (isType<wallFvPatch>(currPatch))
{
const scalarField& nuw = nutb.boundaryField()[patchi];
scalarField magFaceGradU = mag(U.boundaryField()[patchi].snGrad());
forAll(currPatch, facei)
{
label faceCelli = currPatch.faceCells()[facei];
scalar yPlus =
Cmu25*y[patchi][facei]
*::sqrt(k[faceCelli])
/nub.value();
// For corner cells (with two boundary or more faces),
// epsilon and G in the near-wall cell are calculated
// as an average
cellBoundaryFaceCount[faceCelli]++;
epsilon[faceCelli] +=
Cmu75*::pow(k[faceCelli], 1.5)
/(kappa*y[patchi][facei]);
if (yPlus > 11.6)
{
G[faceCelli] +=
nuw[facei]*magFaceGradU[facei]
*Cmu25*::sqrt(k[faceCelli])
/(kappa*y[patchi][facei]);
}
}
}
}
// perform the averaging
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isType<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
}
}
}
}

38
applications/solvers/multiphase/twoPhaseEulerFoam/turbulenceModel/wallViscosity.H
View File

@ -1,38 +0,0 @@
{
scalar Cmu25 = ::pow(Cmu, 0.25);
const fvPatchList& patches = mesh.boundary();
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (isType<wallFvPatch>(currPatch))
{
scalarField& nutw = nutb.boundaryField()[patchi];
forAll(currPatch, facei)
{
label faceCelli = currPatch.faceCells()[facei];
// calculate yPlus
scalar yPlus =
Cmu25*y[patchi][facei]
*::sqrt(k[faceCelli])
/nub.value();
if (yPlus > 11.6)
{
nutw[facei] =
yPlus*nub.value()*kappa
/::log(E*yPlus)
- nub.value();
}
else
{
nutw[facei] = 0.0;
}
}
}
}
}