twoPhaseEulerFoam: Files transferred from bubbleFoam

This commit is contained in:
Henry
2012-11-14 15:45:19 +00:00
parent 63a406e935
commit 8f45806c1d
5 changed files with 409 additions and 0 deletions

View File

@ -0,0 +1,178 @@
IOdictionary RASProperties
(
IOobject
(
"RASProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
Switch turbulence
(
RASProperties.lookup("turbulence")
);
dictionary kEpsilonDict
(
RASProperties.subDictPtr("kEpsilonCoeffs")
);
dimensionedScalar Cmu
(
dimensionedScalar::lookupOrAddToDict
(
"Cmu",
kEpsilonDict,
0.09
)
);
dimensionedScalar C1
(
dimensionedScalar::lookupOrAddToDict
(
"C1",
kEpsilonDict,
1.44
)
);
dimensionedScalar C2
(
dimensionedScalar::lookupOrAddToDict
(
"C2",
kEpsilonDict,
1.92
)
);
dimensionedScalar alpha1k
(
dimensionedScalar::lookupOrAddToDict
(
"alpha1k",
kEpsilonDict,
1.0
)
);
dimensionedScalar alpha1Eps
(
dimensionedScalar::lookupOrAddToDict
(
"alpha1Eps",
kEpsilonDict,
0.76923
)
);
dictionary wallFunctionDict
(
RASProperties.subDictPtr("wallFunctionCoeffs")
);
dimensionedScalar kappa
(
dimensionedScalar::lookupOrAddToDict
(
"kappa",
wallFunctionDict,
0.41
)
);
dimensionedScalar E
(
dimensionedScalar::lookupOrAddToDict
(
"E",
wallFunctionDict,
9.8
)
);
if (RASProperties.lookupOrDefault("printCoeffs", false))
{
Info<< "kEpsilonCoeffs" << kEpsilonDict << nl
<< "wallFunctionCoeffs" << wallFunctionDict << endl;
}
nearWallDist y(mesh);
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
);
Info<< "Calculating field nut2\n" << endl;
volScalarField nut2
(
IOobject
(
"nut2",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Cmu*sqr(k)/epsilon
);
Info<< "Calculating field nuEff1\n" << endl;
volScalarField nuEff1
(
IOobject
(
"nuEff1",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
sqr(Ct)*nut2 + nu1
);
Info<< "Calculating field nuEff2\n" << endl;
volScalarField nuEff2
(
IOobject
(
"nuEff2",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
nut2 + nu2
);

View File

@ -0,0 +1,64 @@
if (turbulence)
{
if (mesh.changing())
{
y.correct();
}
tmp<volTensorField> tgradU2 = fvc::grad(U2);
volScalarField G(2*nut2*(tgradU2() && dev(symm(tgradU2()))));
tgradU2.clear();
#include "wallFunctions.H"
// Dissipation equation
fvScalarMatrix epsEqn
(
fvm::ddt(epsilon)
+ fvm::div(phi2, epsilon)
- fvm::Sp(fvc::div(phi2), epsilon)
- fvm::laplacian
(
alpha1Eps*nuEff2, epsilon,
"laplacian(DepsilonEff,epsilon)"
)
==
C1*G*epsilon/k
- fvm::Sp(C2*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(k)
+ fvm::div(phi2, k)
- fvm::Sp(fvc::div(phi2), k)
- fvm::laplacian
(
alpha1k*nuEff2, k,
"laplacian(DkEff,k)"
)
==
G
- fvm::Sp(epsilon/k, k)
);
kEqn.relax();
kEqn.solve();
k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8));
//- Re-calculate turbulence viscosity
nut2 = Cmu*sqr(k)/epsilon;
#include "wallViscosity.H"
}
nuEff2 = nut2 + nu2;

View File

@ -0,0 +1,50 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 3 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, see <http://www.gnu.org/licenses/>.
Global
wallDissipation
Description
Set wall dissipation in the epsilon matrix
\*---------------------------------------------------------------------------*/
{
const fvPatchList& patches = mesh.boundary();
forAll(patches, patchi)
{
const fvPatch& p = patches[patchi];
if (isA<wallFvPatch>(p))
{
epsEqn.setValues
(
p.faceCells(),
epsilon.boundaryField()[patchi].patchInternalField()
);
}
}
}
// ************************************************************************* //

View File

@ -0,0 +1,81 @@
{
labelList cellBoundaryFaceCount(epsilon.size(), 0);
scalar Cmu25 = ::pow(Cmu.value(), 0.25);
scalar Cmu75 = ::pow(Cmu.value(), 0.75);
scalar kappa_ = kappa.value();
scalar nu2_ = nu2.value();
const fvPatchList& patches = mesh.boundary();
//- Initialise the near-wall P field to zero
forAll(patches, patchi)
{
const fvPatch& currPatch = patches[patchi];
if (isA<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 (isA<wallFvPatch>(currPatch))
{
const scalarField& nut2w = nut2.boundaryField()[patchi];
scalarField magFaceGradU(mag(U2.boundaryField()[patchi].snGrad()));
forAll(currPatch, facei)
{
label faceCelli = currPatch.faceCells()[facei];
// 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]);
G[faceCelli] +=
(nut2w[facei] + nu2_)*magFaceGradU[facei]
*Cmu25*::sqrt(k[faceCelli])
/(kappa_*y[patchi][facei]);
}
}
}
// perform the averaging
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
cellBoundaryFaceCount[faceCelli] = 1;
}
}
}
}

View File

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