Commit 808ea384 authored by Henry's avatar Henry
Browse files

TurbulenceModels/phaseIncompressible/RAS/mixtureKEpsilon/mixtureKEpsilon: preliminary version

parent 0b416be9
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013 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 "mixtureKEpsilon.H"
#include "bound.H"
#include "twoPhaseSystem.H"
#include "fixedValueFvPatchFields.H"
#include "inletOutletFvPatchFields.H"
#include "fvmSup.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RASModels
{
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
mixtureKEpsilon<BasicTurbulenceModel>::mixtureKEpsilon
(
const alphaField& alpha,
const rhoField& rho,
const volVectorField& U,
const surfaceScalarField& alphaPhi,
const surfaceScalarField& phi,
const transportModel& transport,
const word& propertiesName,
const word& type
)
:
eddyViscosity<RASModel<BasicTurbulenceModel> >
(
type,
alpha,
rho,
U,
alphaPhi,
phi,
transport,
propertiesName
),
liquidTurbulencePtr_(NULL),
Cmu_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cmu",
this->coeffDict_,
0.09
)
),
C1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C1",
this->coeffDict_,
1.44
)
),
C2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C2",
this->coeffDict_,
1.92
)
),
C3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C3",
this->coeffDict_,
C2_.value()
)
),
Cp_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cp",
this->coeffDict_,
0.25
)
),
sigmak_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmak",
this->coeffDict_,
1.0
)
),
sigmaEps_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmaEps",
this->coeffDict_,
1.3
)
),
k_
(
IOobject
(
IOobject::groupName("k", U.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
),
epsilon_
(
IOobject
(
IOobject::groupName("epsilon", U.group()),
this->runTime_.timeName(),
this->mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
this->mesh_
)
{
bound(k_, this->kMin_);
bound(epsilon_, this->epsilonMin_);
if (type == typeName)
{
this->printCoeffs(type);
}
}
template<class BasicTurbulenceModel>
wordList mixtureKEpsilon<BasicTurbulenceModel>::epsilonBoundaryTypes
(
const volScalarField& epsilon
) const
{
const volScalarField::GeometricBoundaryField& ebf = epsilon.boundaryField();
wordList ebt = ebf.types();
forAll(ebf, patchi)
{
if (isA<fixedValueFvPatchScalarField>(ebf[patchi]))
{
ebt[patchi] = fixedValueFvPatchScalarField::typeName;
}
}
return ebt;
}
template<class BasicTurbulenceModel>
void mixtureKEpsilon<BasicTurbulenceModel>::correctInletOutlet
(
volScalarField& vsf,
const volScalarField& refVsf
) const
{
volScalarField::GeometricBoundaryField& bf = vsf.boundaryField();
const volScalarField::GeometricBoundaryField& refBf =
refVsf.boundaryField();
forAll(bf, patchi)
{
if
(
isA<inletOutletFvPatchScalarField>(bf[patchi])
&& isA<inletOutletFvPatchScalarField>(refBf[patchi])
)
{
refCast<inletOutletFvPatchScalarField>
(bf[patchi]).refValue() =
refCast<const inletOutletFvPatchScalarField>
(refBf[patchi]).refValue();
}
}
}
template<class BasicTurbulenceModel>
void mixtureKEpsilon<BasicTurbulenceModel>::initMixtureFields()
{
if (rhom_.valid()) return;
// Local references to gas-phase properties
const volScalarField& kg = this->k_;
const volScalarField& epsilong = this->epsilon_;
// Local references to liquid-phase properties
mixtureKEpsilon<BasicTurbulenceModel>& turbc = this->liquidTurbulence();
const volScalarField& kl = turbc.k_;
const volScalarField& epsilonl = turbc.epsilon_;
Ct2_.set
(
new volScalarField
(
IOobject
(
"Ct2",
this->runTime_.timeName(),
this->mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
Ct2()
)
);
rhom_.set
(
new volScalarField
(
IOobject
(
"rhom",
this->runTime_.timeName(),
this->mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
rhom()
)
);
km_.set
(
new volScalarField
(
IOobject
(
"km",
this->runTime_.timeName(),
this->mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mix(kl, kg),
kl.boundaryField().types()
)
);
correctInletOutlet(km_(), kl);
epsilonm_.set
(
new volScalarField
(
IOobject
(
"epsilonm",
this->runTime_.timeName(),
this->mesh_,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mix(epsilonl, epsilong),
epsilonBoundaryTypes(epsilonl)
)
);
correctInletOutlet(epsilonm_(), epsilonl);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
bool mixtureKEpsilon<BasicTurbulenceModel>::read()
{
if (eddyViscosity<RASModel<BasicTurbulenceModel> >::read())
{
Cmu_.readIfPresent(this->coeffDict());
C1_.readIfPresent(this->coeffDict());
C2_.readIfPresent(this->coeffDict());
C3_.readIfPresent(this->coeffDict());
Cp_.readIfPresent(this->coeffDict());
sigmak_.readIfPresent(this->coeffDict());
sigmaEps_.readIfPresent(this->coeffDict());
return true;
}
else
{
return false;
}
}
template<class BasicTurbulenceModel>
void mixtureKEpsilon<BasicTurbulenceModel>::correctNut()
{
this->nut_ = Cmu_*sqr(k_)/epsilon_;
this->nut_.correctBoundaryConditions();
}
template<class BasicTurbulenceModel>
mixtureKEpsilon<BasicTurbulenceModel>&
mixtureKEpsilon<BasicTurbulenceModel>::liquidTurbulence() const
{
if (!liquidTurbulencePtr_)
{
const volVectorField& U = this->U_;
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
const transportModel& liquid = fluid.otherPhase(gas);
liquidTurbulencePtr_ =
&const_cast<mixtureKEpsilon<BasicTurbulenceModel>&>
(
U.db().lookupObject<mixtureKEpsilon<BasicTurbulenceModel> >
(
IOobject::groupName
(
turbulenceModel::propertiesName,
liquid.name()
)
)
);
}
return *liquidTurbulencePtr_;
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::Ct2() const
{
const mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =
this->liquidTurbulence();
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
const transportModel& liquid = fluid.otherPhase(gas);
const volScalarField& alphag = this->alpha_;
volScalarField magUr(mag(liquidTurbulence.U() - this->U()));
volScalarField beta
(
(6*this->Cmu_/(4*sqrt(3.0/2.0)))
*alphag*fluid.drag(gas).K(magUr)/liquid.rho()
*(liquidTurbulence.k_/liquidTurbulence.epsilon_)
);
volScalarField Ct0((3 + beta)/(1 + beta + 2*gas.rho()/liquid.rho()));
volScalarField fAlphad((180 + (-4.71e3 + 4.26e4*alphag)*alphag)*alphag);
return sqr(1 + (Ct0 - 1)*exp(-fAlphad));
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rholEff() const
{
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
return fluid.otherPhase(gas).rho();
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rhogEff() const
{
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
return
gas.rho()
+ fluid.Cvm()*fluid.otherPhase(gas).rho();
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::rhom() const
{
const volScalarField& alphag = this->alpha_;
const volScalarField& alphal = this->liquidTurbulence().alpha_;
return alphal*rholEff() + alphag*rhogEff();
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mix
(
const volScalarField& fc,
const volScalarField& fd
) const
{
const volScalarField& alphag = this->alpha_;
const volScalarField& alphal = this->liquidTurbulence().alpha_;
return (alphal*rholEff()*fc + alphag*rhogEff()*fd)/rhom_();
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mixU
(
const volScalarField& fc,
const volScalarField& fd
) const
{
const volScalarField& alphag = this->alpha_;
const volScalarField& alphal = this->liquidTurbulence().alpha_;
return
(alphal*rholEff()*fc + alphag*rhogEff()*Ct2_()*fd)
/(alphal*rholEff() + alphag*rhogEff()*Ct2_());
}
template<class BasicTurbulenceModel>
tmp<surfaceScalarField> mixtureKEpsilon<BasicTurbulenceModel>::mixFlux
(
const surfaceScalarField& fc,
const surfaceScalarField& fd
) const
{
const volScalarField& alphag = this->alpha_;
const volScalarField& alphal = this->liquidTurbulence().alpha_;
surfaceScalarField alphalf(fvc::interpolate(alphal));
surfaceScalarField alphagf(fvc::interpolate(alphag));
surfaceScalarField rholEfff(fvc::interpolate(rholEff()));
surfaceScalarField rhogEfff(fvc::interpolate(rhogEff()));
return
fvc::interpolate(rhom_()/(alphal*rholEff() + alphag*rhogEff()*Ct2_()))
*(alphalf*rholEfff*fc + alphagf*rhogEfff*fvc::interpolate(Ct2_())*fd);
}
template<class BasicTurbulenceModel>
tmp<volScalarField> mixtureKEpsilon<BasicTurbulenceModel>::bubbleG() const
{
const mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =
this->liquidTurbulence();
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
const transportModel& liquid = fluid.otherPhase(gas);
volScalarField magUr(mag(liquidTurbulence.U() - this->U()));
// Lahey model
tmp<volScalarField> bubbleG
(
Cp_
*sqr(liquid)*liquid.rho()
*(
pow3(magUr)
+ pow(fluid.drag(gas).K(magUr)*gas.d()/liquid.rho(), 3.0/4.0)
*pow(magUr, 9.0/4.0)
)
*gas
/gas.d()
);
// Simple model
// tmp<volScalarField> bubbleG
// (
// Cp_*sqr(liquid)*gas*fluid.drag(gas).K(magUr)*sqr(magUr)
// );
return bubbleG;
}
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::kSource() const
{
return fvm::Su(bubbleG(), km_());
}
template<class BasicTurbulenceModel>
tmp<fvScalarMatrix> mixtureKEpsilon<BasicTurbulenceModel>::epsilonSource() const
{
return fvm::Su(C3_*epsilonm_()*bubbleG()/km_(), epsilonm_());
}
template<class BasicTurbulenceModel>
void mixtureKEpsilon<BasicTurbulenceModel>::correct()
{
const transportModel& gas = this->transport();
const twoPhaseSystem& fluid = gas.fluid();
// Only solve the mixture turbulence for the gas-phase
if (&gas != &fluid.phase1())
{
// This is the liquid phase but check the model for the gas-phase
// is consistent
this->liquidTurbulence();
return;
}
if (!this->turbulence_)
{
return;
}
// Initialise the mixture fields if they have not yet been constructed
initMixtureFields();
// Local references to gas-phase properties
const surfaceScalarField& phig = this->phi_;
const volVectorField& Ug = this->U_;
const volScalarField& alphag = this->alpha_;
volScalarField& kg = this->k_;
volScalarField& epsilong = this->epsilon_;
volScalarField& nutg = this->nut_;
// Local references to liquid-phase properties
mixtureKEpsilon<BasicTurbulenceModel>& liquidTurbulence =
this->liquidTurbulence();
const surfaceScalarField& phil = liquidTurbulence.phi_;
const volVectorField& Ul = liquidTurbulence.U_;
const volScalarField& alphal = liquidTurbulence.alpha_;
volScalarField& kl = liquidTurbulence.k_;
volScalarField& epsilonl = liquidTurbulence.epsilon_;
volScalarField& nutl = liquidTurbulence.nut_;
// Local references to mixture properties
volScalarField& rhom = rhom_();
volScalarField& km = km_();
volScalarField& epsilonm = epsilonm_();
eddyViscosity<RASModel<BasicTurbulenceModel> >::correct();
// Update the effective mixture density
rhom = this->rhom();
// Mixture flux
surfaceScalarField phim("phim", mixFlux(phil, phig));
// Mixture velocity divergence
volScalarField divUm
(
mixU
(
fvc::div(fvc::absolute(phil, Ul)),
fvc::div(fvc::absolute(phig, Ug))
)
);
tmp<volScalarField> Gc;
{
tmp<volTensorField> tgradUl =