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

driftFluxFoam: Change turbulence modelling to new templated framework

Browse Source 
Added k-epsilon model with buoyancy correction
This commit is contained in:
Henry
2014-04-28 15:31:28 +01:00
committed by Andrew Heather
parent 6723be5056
commit 63bef62f8e
33 changed files with 668 additions and 618 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/solvers/multiphase/driftFluxFoam/Make/files
View File

@ -1,4 +1,5 @@
incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.C
compressibleTurbulenceModels.C
driftFluxFoam.C
EXE = $(FOAM_APPBIN)/driftFluxFoam

8
applications/solvers/multiphase/driftFluxFoam/Make/options
View File

@ -9,13 +9,17 @@ EXE_INC = \
-I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I./relativeVelocityModels/lnInclude
EXE_LIBS = \
-ldriftFluxTransportModels \
-ldriftFluxRelativeVelocityModels \
-lfiniteVolume \
-lmeshTools \
-lsampling \
-lfvOptions \
-lincompressibleTransportModels \
-ldriftFluxTransportModels \
-ldriftFluxRelativeVelocityModels
-lturbulenceModels \
-lcompressibleTurbulenceModels

3
applications/solvers/multiphase/driftFluxFoam/UEqn.H
View File

@ -5,8 +5,7 @@
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
+ fvc::div(UdmModel.tauDm())
- fvm::laplacian(muEff, U)
- fvc::div(muEff*dev2(T(fvc::grad(U))))
+ turbulence->divDevRhoReff(U)
==
fvOptions(rho, U)
);

2
applications/solvers/multiphase/driftFluxFoam/alphaEqn.H
View File

@ -27,7 +27,7 @@
mesh,
linear<scalar>(mesh),
fv::uncorrectedSnGrad<scalar>(mesh)
).fvmLaplacian(fvc::interpolate(mut/rho), alpha1)
).fvmLaplacian(fvc::interpolate(turbulence->nut()), alpha1)
);
Info<< "Phase-1 volume fraction = "

4
applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H
View File

@ -55,7 +55,7 @@
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1) - fvc::ddt(alpha1)
- fvm::laplacian(mut/rho, alpha1)
- fvm::laplacian(turbulence->nut(), alpha1)
);
alpha1Eqn.solve(mesh.solver("alpha1Diffusion"));
@ -71,5 +71,5 @@
}
rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
rho == alpha1*rho1 + alpha2*rho2;
rho = mixture.rho();
}

73
applications/solvers/multiphase/driftFluxFoam/compressibleTurbulenceModels.C Normal file
View File

@ -0,0 +1,73 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 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/>.
\*---------------------------------------------------------------------------*/
#include "CompressibleTurbulenceModel.H"
#include "incompressibleTwoPhaseInteractingMixture.H"
#include "addToRunTimeSelectionTable.H"
#include "makeTurbulenceModel.H"
#include "laminar.H"
#include "RASModel.H"
#include "LESModel.H"
makeBaseTurbulenceModel
(
geometricOneField,
volScalarField,
compressibleTurbulenceModel,
CompressibleTurbulenceModel,
incompressibleTwoPhaseInteractingMixture
);
#define makeRASModel(Type) \
makeTemplatedTurbulenceModel \
( \
incompressibleTwoPhaseInteractingMixtureCompressibleTurbulenceModel, \
RAS, \
Type \
)
#define makeLESModel(Type) \
makeTemplatedTurbulenceModel \
( \
incompressibleTwoPhaseInteractingMixtureCompressibleTurbulenceModel, \
LES, \
Type \
)
#include "kEpsilon.H"
makeRASModel(kEpsilon);
#include "buoyantKEpsilon.H"
makeRASModel(buoyantKEpsilon);
#include "Smagorinsky.H"
makeLESModel(Smagorinsky);
#include "kEqn.H"
makeLESModel(kEqn);
// ************************************************************************* //

196
applications/solvers/multiphase/driftFluxFoam/createFields.H
View File

@ -29,17 +29,14 @@
#include "createPhi.H"
// Transport
// ~~~~~~~~~
Info<< "Reading transportProperties\n" << endl;
incompressibleTwoPhaseInteractingMixture twoPhaseProperties(U, phi);
incompressibleTwoPhaseInteractingMixture mixture(U, phi);
volScalarField& alpha1(twoPhaseProperties.alpha1());
volScalarField& alpha2(twoPhaseProperties.alpha2());
volScalarField& alpha1(mixture.alpha1());
volScalarField& alpha2(mixture.alpha2());
const dimensionedScalar& rho1 = twoPhaseProperties.rhod();
const dimensionedScalar& rho2 = twoPhaseProperties.rhoc();
const dimensionedScalar& rho1 = mixture.rhod();
const dimensionedScalar& rho2 = mixture.rhoc();
IOdictionary transportProperties
(
@ -64,9 +61,8 @@
IOobject::NO_READ,
IOobject::NO_WRITE
),
alpha1*rho1 + alpha2*rho2
mixture.rho()
);
rho.oldTime();
// Mass flux
surfaceScalarField rhoPhi
@ -84,190 +80,28 @@
// Relative Velocity
// ~~~~~~~~~~~~~~~~~
autoPtr<relativeVelocityModel> UdmModelPtr
(
relativeVelocityModel::New
(
transportProperties,
twoPhaseProperties
mixture
)
);
relativeVelocityModel& UdmModel(UdmModelPtr());
// Turbulence
// ~~~~~~~~~~
IOdictionary RASProperties
// Construct compressible turbulence model
autoPtr
<
CompressibleTurbulenceModel<incompressibleTwoPhaseInteractingMixture>
> turbulence
(
IOobject
(
"RASProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
CompressibleTurbulenceModel<incompressibleTwoPhaseInteractingMixture>
::New(rho, U, rhoPhi, mixture)
);
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 C3
(
dimensionedScalar::lookupOrAddToDict
(
"C3",
kEpsilonDict,
0.85
)
);
dimensionedScalar sigmak
(
dimensionedScalar::lookupOrAddToDict
(
"sigmak",
kEpsilonDict,
1.0
)
);
dimensionedScalar sigmaEps
(
dimensionedScalar::lookupOrAddToDict
(
"sigmaEps",
kEpsilonDict,
1.3
)
);
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 mut\n" << endl;
volScalarField mut
(
IOobject
(
"mut",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Cmu*rho*sqr(k)/epsilon
);
Info<< "Calculating field muEff\n" << endl;
volScalarField muEff
(
IOobject
(
"muEff",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mut + twoPhaseProperties.mu()
);
// Pressure
// ~~~~~~~~
Info<< "Calculating field (g.h)f\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
@ -309,6 +143,4 @@
// MULES Correction
// ~~~~~~~~~~~~~~~~
tmp<surfaceScalarField> tphiAlphaCorr0;

10
applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C
View File

@ -38,10 +38,8 @@ Description
#include "subCycle.H"
#include "incompressibleTwoPhaseInteractingMixture.H"
#include "relativeVelocityModel.H"
#include "nearWallDist.H"
#include "wallFvPatch.H"
#include "bound.H"
#include "Switch.H"
#include "turbulenceModel.H"
#include "CompressibleTurbulenceModel.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "fixedFluxPressureFvPatchScalarField.H"
@ -86,7 +84,7 @@ int main(int argc, char *argv[])
#include "alphaEqnSubCycle.H"
twoPhaseProperties.correct();
mixture.correct();
#include "UEqn.H"
@ -98,7 +96,7 @@ int main(int argc, char *argv[])
if (pimple.turbCorr())
{
#include "kEpsilon.H"
turbulence->correct();
}
}

29
applications/solvers/multiphase/driftFluxFoam/incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.H
View File

@ -133,21 +133,36 @@ public:
return mu_;
}
//- Return the dynamic mixture viscosity for patch
virtual tmp<scalarField> mu(const label patchi) const
{
return mu_.boundaryField()[patchi];
}
//- Return the mixture density
virtual tmp<volScalarField> rho() const
{
return alpha1_*rhod_ + alpha2_*rhoc_;
}
//- Return the mixture density for patch
virtual tmp<scalarField> rho(const label patchi) const
{
return
alpha1_.boundaryField()[patchi]*rhod_.value()
+ alpha2_.boundaryField()[patchi]*rhoc_.value();
}
//- Return the mixture viscosity
virtual tmp<volScalarField> nu() const
{
notImplemented("incompressibleTwoPhaseInteractingMixture::nu()");
return volScalarField::null();
return mu_/rho();
}
//- Return the mixture viscosity for patch
virtual tmp<scalarField> nu(const label patchi) const
{
notImplemented
(
"incompressibleTwoPhaseInteractingMixture::nu(const label)"
);
return scalarField::null();
return mu_.boundaryField()[patchi]/rho(patchi);
}
//- Correct the laminar viscosity

82
applications/solvers/multiphase/driftFluxFoam/kEpsilon.H
View File

@ -1,82 +0,0 @@
if (turbulence)
{
if (mesh.changing())
{
y.correct();
}
dimensionedScalar k0("k0", k.dimensions(), 0);
dimensionedScalar kMin("kMin", k.dimensions(), SMALL);
dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0);
dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField G(mut*(tgradU() && dev(twoSymm(tgradU()))));
tgradU.clear();
volScalarField Gcoef
(
Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin)
);
volScalarField mul(twoPhaseProperties.mu());
#include "wallFunctions.H"
// Dissipation equation
fvScalarMatrix epsEqn
(
fvm::ddt(rho, epsilon)
+ fvm::div(rhoPhi, epsilon)
- fvm::laplacian
(
mut/sigmaEps + mul, epsilon,
"laplacian(DepsilonEff,epsilon)"
)
==
C1*G*epsilon/(k + kMin)
- fvm::SuSp(C1*(1.0 - C3)*Gcoef, epsilon)
- fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon)
);
#include "wallDissipation.H"
epsEqn.relax();
epsEqn.solve();
bound(epsilon, epsilon0);
// Turbulent kinetic energy equation
fvScalarMatrix kEqn
(
fvm::ddt(rho, k)
+ fvm::div(rhoPhi, k)
- fvm::laplacian
(
mut/sigmak + mul, k,
"laplacian(DkEff,k)"
)
==
G
- fvm::SuSp(Gcoef, k)
- fvm::Sp(rho*epsilon/(k + kMin), k)
);
kEqn.relax();
kEqn.solve();
bound(k, k0);
//- Re-calculate viscosity
mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin);
#include "wallViscosity.H"
muEff = mut + mul;
}
else
{
muEff = mut + twoPhaseProperties.mu();
}

50
applications/solvers/multiphase/driftFluxFoam/wallDissipation.H
View File

@ -1,50 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 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()
);
}
}
}
// ************************************************************************* //

85
applications/solvers/multiphase/driftFluxFoam/wallFunctions.H
View File

@ -1,85 +0,0 @@
{
labelList cellBoundaryFaceCount(epsilon.size(), 0);
const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
const scalar Cmu75 = ::pow(Cmu.value(), 0.75);
const scalar kappa_ = kappa.value();
const fvPatchList& patches = mesh.boundary();
//- Initialise the near-wall P field to zero
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
forAll(curPatch, facei)
{
label faceCelli = curPatch.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& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
const scalarField& mutw = mut.boundaryField()[patchi];
const scalarField& mulw = mul.boundaryField()[patchi];
scalarField magFaceGradU
(
mag(U.boundaryField()[patchi].snGrad())
);
forAll(curPatch, facei)
{
label faceCelli = curPatch.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] +=
(mutw[facei] + mulw[facei])
*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;
}
}
}
}

38
applications/solvers/multiphase/driftFluxFoam/wallViscosity.H
View File

@ -1,38 +0,0 @@
{
const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
const scalar kappa_ = kappa.value();
const scalar E_ = E.value();
const fvPatchList& patches = mesh.boundary();
forAll(patches, patchi)
{
const fvPatch& curPatch = patches[patchi];
if (isA<wallFvPatch>(curPatch))
{
scalarField& mutw = mut.boundaryField()[patchi];
const scalarField& mulw = mul.boundaryField()[patchi];
const scalarField& rhow = rho.boundaryField()[patchi];
forAll(curPatch, facei)
{
label faceCelli = curPatch.faceCells()[facei];
scalar yPlus =
Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
/(mulw[facei]/rhow[facei]);
if (yPlus > 11.6)
{
mutw[facei] =
mulw[facei]*(yPlus*kappa_/::log(E_*yPlus) - 1);
}
else
{
mutw[facei] = 0.0;
}
}
}
}
}