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 of Foam::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.

Admin message