[Hamed1117]

Hello everyone,

Currently I am working on implementation of a non-Newtonian viscosity model in src/thermophysicalModel which requires some code modifications in src/thermophysicalMode/basic and src/thermophysicalMode/specie and I am using the following paper (https://doi.org/10.4236/ojfd.2020.102009) as the initial reference.

The new viscosity model needs to be modified inside specie/newNon-Newtonian.C and called by basic/rhoThermo/heRhoThermo.C where the shear rate (sr) and phase volume fraction (alpha.1) are the main variables.

The new model is compiled without any major problem, but some problems arises when it comes to the calling part using heRhoThermo.C.

According to the paper, I need to use the objectRegistery tool in order to get access to U1(for sr_) and alpha.1 fields.



heRhoThermo.C

#include "heRhoThermo.H"
#include "fvc.H"
//#include "fvCFD.H"
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

template<class BasicPsiThermo, class MixtureType>
void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::calculate()
{

    const scalarField& hCells = this->he();
    const scalarField& pCells = this->p_;

    scalarField& TCells = this->T_.primitiveFieldRef();
    scalarField& psiCells = this->psi_.primitiveFieldRef();
    scalarField& rhoCells = this->rho_.primitiveFieldRef();
    scalarField& muCells = this->mu_.primitiveFieldRef();
     scalarField& alphaCells = this->alpha_.primitiveFieldRef();


//register alpha.1_
    const volScalarField& alpha1 =  this->db().objectRegistry::lookupObject<volScalarField>("alpha.particles");
    volScalarField alpha1_ = alpha1;
    scalarField& alpha1Cells = alpha1_.primitiveFieldRef();

//register U1 for sr_ 
 const volVectorField& U1 = this->db().objectRegistry::lookupObject<volVectorField>("U.particles");    
 volScalarField sr_=sqrt(2.0)*mag(symm(fvc::grad(U1)));
  scalarField& srCells = sr_.primitiveFieldRef();
   
 forAll(TCells, celli)
    {
        const typename MixtureType::thermoType& mixture_ =
            this->cellMixture(celli);


        TCells[celli] = mixture_.THE
        (
            hCells[celli],
            pCells[celli],
            TCells[celli]
        );

        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
        rhoCells[celli] = mixture_.rho(pCells[celli], TCells[celli]);
        muCells[celli] = mixture_.mu(pCells[celli], TCells[celli], alpha1Cells[celli], srCells[celli]);
        alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli], alpha1Cells[celli], srCells[celli]);
    }

    volScalarField::Boundary& pBf =
        this->p_.boundaryFieldRef();

    volScalarField::Boundary& TBf =
        this->T_.boundaryFieldRef();
//
    volScalarField::Boundary& alpha1Bf =
             alpha1_.boundaryFieldRef();

    volScalarField::Boundary& srBf =
             sr_.boundaryFieldRef();
//
    volScalarField::Boundary& psiBf =
        this->psi_.boundaryFieldRef();
volScalarField::Boundary& psiBf =
        this->psi_.boundaryFieldRef();

    volScalarField::Boundary& rhoBf =
        this->rho_.boundaryFieldRef();

    volScalarField::Boundary& heBf =
        this->he().boundaryFieldRef();

    volScalarField::Boundary& muBf =
        this->mu_.boundaryFieldRef();

    volScalarField::Boundary& alphaBf =
        this->alpha_.boundaryFieldRef();

    forAll(this->T_.boundaryField(), patchi)
    {
        fvPatchScalarField& pp = pBf[patchi];
        fvPatchScalarField& pT = TBf[patchi];
//
        fvPatchScalarField& palpha1 = alpha1Bf[patchi];
        fvPatchScalarField& psr = srBf[patchi];
//
        fvPatchScalarField& ppsi = psiBf[patchi];
        fvPatchScalarField& prho = rhoBf[patchi];
        fvPatchScalarField& phe = heBf[patchi];
        fvPatchScalarField& pmu = muBf[patchi];
        fvPatchScalarField& palpha = alphaBf[patchi];
        if (pT.fixesValue())
        {

            forAll(pT, facei)
            {
                const typename MixtureType::thermoType& mixture_ =
                    this->patchFaceMixture(patchi, facei);
               phe[facei] = mixture_.HE(pp[facei], pT[facei]);
                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
                prho[facei] = mixture_.rho(pp[facei], pT[facei]);
                pmu[facei] = mixture_.mu(pp[facei], pT[facei],palpha1[facei], psr[facei]);
                palpha[facei] = mixture_.alphah(pp[facei], pT[facei], palpha1[facei], psr[facei]);
            }
        }
        else
        {
            forAll(pT, facei)
            {
                const typename MixtureType::thermoType& mixture_ =
                    this->patchFaceMixture(patchi, facei);
                pT[facei] = mixture_.THE(phe[facei], pp[facei], pT[facei]);
                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
                prho[facei] = mixture_.rho(pp[facei], pT[facei]);
//
                pmu[facei] = mixture_.mu(pp[facei], pT[facei], palpha1[facei], psr[facei]);
//
                palpha[facei] = mixture_.alphah(pp[facei], pT[facei], palpha1[facei], psr[facei]);
            }
        }
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template<class BasicPsiThermo, class MixtureType>
Foam::heRhoThermo<BasicPsiThermo, MixtureType>::heRhoThermo
(
    const fvMesh& mesh,
    const word& phaseName
)
:
    heThermo<BasicPsiThermo, MixtureType>(mesh, phaseName)
{
    calculate();
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

template<class BasicPsiThermo, class MixtureType>
Foam::heRhoThermo<BasicPsiThermo, MixtureType>::~heRhoThermo()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template<class BasicPsiThermo, class MixtureType>
void Foam::heRhoThermo<BasicPsiThermo, MixtureType>::correct()
{
    if (debug)
    {
        InfoInFunction << endl;
    }

    calculate();

    if (debug)
    {
        Info<< "    Finished" << endl;
    }
}

And it compiles successfully. But when I run the solver in a simple case I got the following error as:

--> FOAM FATAL ERROR: 

    request for volVectorField U.particles from objectRegistry region0 failed
    available objects of type volVectorField are
0()

    From function const Type& Foam::objectRegistry::lookupObject(const Foam::word&) const [with Type = Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>]
    in file /opt/openfoam5/src/OpenFOAM/lnInclude/objectRegistryTemplates.C at line 193.

FOAM aborting

It is clear that I make a mistake during the U.particles registration process. The paper says that U (volVectorField) needs to be called in viscosity claculation loop (paper - listing.1), but I wonder where I should add it !! (inside basic/heRothermo.C, mySolver.C or specie/newNonNewtonianModel.C ?? )