6c0087d005
The pressure work term for total internal energy is div(U p) which can be
discretised is various ways, given a mass flux field phi it seems logical to
implement it in the form div(phi/interpolate(rho), p) but this is not exactly
consistent with the relationship between enthalpy and internal energy (h = e +
p/rho) and the transport of enthalpy, it would be more consistent to implement
it in the form div(phi, p/rho). A further improvement in consistency can be
gained by using the same convection scheme for this work term and the convection
term div(phi, e) and for reacting solvers this is easily achieved by using the
multi-variate limiter mvConvection provided for energy and specie convection.
This more consistent total internal energy work term has now been implemented in
all the compressible and reacting flow solvers and provides more accurate
solutions when running with internal energy, particularly for variable density
mixing cases with small pressure variation.
For non-reacting compressible solvers this improvement requires a change to the
corresponding divScheme in fvSchemes:
"div\(alphaPhi.*,p\)" -> "div\(alphaRhoPhi.*,\(p\|thermo:rho.*\)\)"
and all the tutorials have been updated accordingly.
157 lines
4.3 KiB
C++
157 lines
4.3 KiB
C++
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration | Website: https://openfoam.org
|
|
\\ / A nd | Copyright (C) 2015-2021 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 "AnisothermalPhaseModel.H"
|
|
#include "phaseSystem.H"
|
|
|
|
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
|
|
|
|
template<class BasePhaseModel>
|
|
Foam::tmp<Foam::volScalarField>
|
|
Foam::AnisothermalPhaseModel<BasePhaseModel>::filterPressureWork
|
|
(
|
|
const tmp<volScalarField>& pressureWork
|
|
) const
|
|
{
|
|
const volScalarField& alpha = *this;
|
|
|
|
scalar pressureWorkAlphaLimit =
|
|
this->thermo_->properties()
|
|
.lookupOrDefault("pressureWorkAlphaLimit", 0.0);
|
|
|
|
if (pressureWorkAlphaLimit > 0)
|
|
{
|
|
return
|
|
(
|
|
max(alpha - pressureWorkAlphaLimit, scalar(0))
|
|
/max(alpha - pressureWorkAlphaLimit, pressureWorkAlphaLimit)
|
|
)*pressureWork;
|
|
}
|
|
else
|
|
{
|
|
return pressureWork;
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
template<class BasePhaseModel>
|
|
Foam::AnisothermalPhaseModel<BasePhaseModel>::AnisothermalPhaseModel
|
|
(
|
|
const phaseSystem& fluid,
|
|
const word& phaseName,
|
|
const bool referencePhase,
|
|
const label index
|
|
)
|
|
:
|
|
BasePhaseModel(fluid, phaseName, referencePhase, index)
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
|
|
|
template<class BasePhaseModel>
|
|
Foam::AnisothermalPhaseModel<BasePhaseModel>::~AnisothermalPhaseModel()
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
template<class BasePhaseModel>
|
|
void Foam::AnisothermalPhaseModel<BasePhaseModel>::correctThermo()
|
|
{
|
|
BasePhaseModel::correctThermo();
|
|
|
|
this->thermo_->correct();
|
|
}
|
|
|
|
|
|
template<class BasePhaseModel>
|
|
bool Foam::AnisothermalPhaseModel<BasePhaseModel>::isothermal() const
|
|
{
|
|
return false;
|
|
}
|
|
|
|
|
|
template<class BasePhaseModel>
|
|
Foam::tmp<Foam::fvScalarMatrix>
|
|
Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
|
|
{
|
|
const volScalarField& alpha = *this;
|
|
const volScalarField& rho = this->rho();
|
|
|
|
const tmp<volVectorField> tU(this->U());
|
|
const volVectorField& U(tU());
|
|
|
|
const tmp<surfaceScalarField> talphaRhoPhi(this->alphaRhoPhi());
|
|
const surfaceScalarField& alphaRhoPhi(talphaRhoPhi());
|
|
|
|
const tmp<volScalarField> tcontErr(this->continuityError());
|
|
const volScalarField& contErr(tcontErr());
|
|
|
|
tmp<volScalarField> tK(this->K());
|
|
const volScalarField& K(tK());
|
|
|
|
volScalarField& he = this->thermo_->he();
|
|
|
|
tmp<fvScalarMatrix> tEEqn
|
|
(
|
|
fvm::ddt(alpha, rho, he)
|
|
+ fvm::div(alphaRhoPhi, he)
|
|
- fvm::Sp(contErr, he)
|
|
|
|
+ fvc::ddt(alpha, rho, K) + fvc::div(alphaRhoPhi, K)
|
|
- contErr*K
|
|
|
|
+ this->divq(he)
|
|
==
|
|
alpha*this->Qdot()
|
|
);
|
|
|
|
// Add the appropriate pressure-work term
|
|
if (he.name() == this->thermo_->phasePropertyName("e"))
|
|
{
|
|
tEEqn.ref() += filterPressureWork
|
|
(
|
|
fvc::div
|
|
(
|
|
fvc::absolute(alphaRhoPhi, alpha, rho, U),
|
|
this->thermo().p()/rho
|
|
)
|
|
+ (fvc::ddt(alpha) - contErr/rho)*this->thermo().p()
|
|
);
|
|
}
|
|
else if (this->thermo_->dpdt())
|
|
{
|
|
tEEqn.ref() -= filterPressureWork(alpha*this->fluid().dpdt());
|
|
}
|
|
|
|
return tEEqn;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|