buoyancyTurbSource not working with interFoam
Summary
The fvOption buoyancyTurbSource does not work with multiphase solvers like interFoam even though, looking at the code, it seems like this was intended.
Steps to reproduce
- Copy a random turbulent interFoam tutorial. E.g.
tutorials/multiphase/interFoam/RAS/damBreak/damBreak/
- Add a fvOptions file with the following content:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1812 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
BuoyancyTurbSource1
{
// Mandatory entries (unmodifiable)
type buoyancyTurbSource;
active yes;
selectionMode all;
// Optional entries (unmodifiable)
beta 3.3e-03;
rho rho;
alphat alphat;
T T;
}
// ************************************************************************* //
- Run, using the
Allrun
script
What is the current bug behaviour?
You will get the following (or similar) error:
--> FOAM FATAL ERROR: (openfoam-2206)
failed lookup of alphat (objectRegistry region0)
available objects of type volScalarField:
15
(
alpha.water_0
interfaceProperties:K
nut
alpha.water
rho
k
p_rgh
nu
gh
nu1
p
rho_0
nu2
alpha.air
epsilon
)
From const Type& Foam::objectRegistry::lookupObject(const Foam::word&, bool) const [with Type = Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>]
in file /opt/OpenFOAM/OpenFOAM-v2206/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 571.
FOAM exiting
What is the expected correct behavior?
At first, it may seem like the above error is the correct behavior and the implemented buoyancy source terms are simply not meant for isothermal solvers. However, taking a closer look at the code (in src/fvOptions/sources/derived/buoyancyTurbSource/buoyancyTurbSource.C
), it appears that the developers intended this function to be called for isothermal multiphase solvers:
void Foam::fv::buoyancyTurbSource::addSup
(
const volScalarField& alpha,
const volScalarField& rho,
fvMatrix<scalar>& eqn,
const label fieldi
)
{
if (fieldi == 1)
{
buoyancyTurbSourceK(alpha, rho, eqn, fieldi);
return;
}
}
This would call the function defined in src/fvOptions/sources/derived/buoyancyTurbSource/buoyancyTurbSourceTemplates.C
, which would not produce the above error.
In fact this would nicely replicate the method described by Brecht Devolder https://doi.org/10.1016/j.coastaleng.2018.04.011
Relevant logs and/or images
The addSup function that (I think) should be called:
void Foam::fv::buoyancyTurbSource::addSup
(
const volScalarField& alpha,
const volScalarField& rho,
fvMatrix<scalar>& eqn,
const label fieldi
)
{
if (fieldi == 1)
{
buoyancyTurbSourceK(alpha, rho, eqn, fieldi);
return;
}
}
which would invoke:
template<class AlphaFieldType, class RhoFieldType>
void Foam::fv::buoyancyTurbSource::buoyancyTurbSourceK
(
const AlphaFieldType& alpha,
const RhoFieldType& rho,
fvMatrix<scalar>& eqn,
const label fieldi
) const
{
const volScalarField& k = eqn.psi();
const dimensionedScalar k0(k.dimensions(), SMALL);
const auto* turbPtr =
mesh_.findObject<turbulenceModel>
(
turbulenceModel::propertiesName
);
const volScalarField& nut = turbPtr->nut();
const dictionary& turbDict = turbPtr->coeffDict();
const scalar Prt
(
turbDict.getCheckOrDefault<scalar>("Prt", 0.85, scalarMinMax::ge(SMALL))
);
// (DTR:Eq. 21, buoyancy correction term)
const tmp<volScalarField> GbByK((nut/Prt)*(fvc::grad(rho) & g_)/max(k, k0));
eqn -= fvm::Sp(GbByK, k);
}
The addSup function that is (apparently) called:
void Foam::fv::buoyancyTurbSource::addSup
(
fvMatrix<scalar>& eqn,
const label fieldi
)
{
if (fieldi == 1)
{
buoyancyTurbSourceK(eqn);
return;
}
if (isEpsilon_)
{
buoyancyTurbSourceEpsilon(eqn);
}
else
{
buoyancyTurbSourceOmega(eqn);
}
}
Environment information
- OpenFOAM version : v2206
- Operating system : openSUSE Leap 15.4
- Hardware info : not relevant
- Compiler : gcc
Possible fixes
I am sure somebody else will have a much clearer idea of how to get the "Devolder version" of buoyancyTurbSource working with interFoam, but here is what I have been able to find out:
- The logic behind the selection of the version of the addSup function to call is found in
src/finiteVolume/cfdTools/general/fvOptions/fvOptionListTemplates.C
and seems quite complicated. To be honest, I do not get which version ofFoam::fv::optionList::operator()
is finally called as the turbulence model adds the source term and which version of addSup is called as a consequence. My feeling is, something is strange there. I have also experienced the source term to have strange units in a different context.
I would be very interested to hear what the issue is (if there is indeed one). It would be very nice to get buoyancyTurbSource working with interFoam, as I had good experiences using Devolder's buoyancy term in the past.