Files
OpenFOAM-12/applications/modules/multiphaseEuler/phaseSystem/phaseModels/AnisothermalPhaseModel/AnisothermalPhaseModel.C
Will Bainbridge 597121a4a7 multiphaseEuler: Library reorganisation
This change makes multiphaseEuler more consistent with other modules and
makes its sub-libraries less inter-dependent. Some left-over references
to multiphaseEulerFoam have also been removed.
2023-09-15 14:45:26 +01:00

163 lines
4.5 KiB
C++

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2015-2023 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"
#include "fvcMeshPhi.H"
#include "fvcDdt.H"
#include "fvmDiv.H"
#include "fvmSup.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),
g_(fluid.mesh().lookupObject<uniformDimensionedVectorField>("g"))
{}
// * * * * * * * * * * * * * * * * 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*rho*(U&g_)
+ 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;
}
// ************************************************************************* //