[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 ?? )

[Hamed1117]

The above mentioned object registery is correct, and it works fine for different solvers (e.g. rhoPimpleFoam, compressibleInterFoam or ...) that are using the same thermoPhysical models. The twoPhaseEulerFoam solver, however, needs a little bit different procedure for accessing the U field.

I described it in the following post:

request fo volVectorField U.particles from objectRegistry regin0 failed-TwoPhaseEuler

[brunomramoa]

Hello,



I am also trying to implement a shear-rate dependent viscosity in the thermophysical models. However, I am getting troubles in the compilation. Would you mind sharing your experience?



Post: Debug help: Viscosity Model in compressibleInterFoam

[hebe03364]

Hello!


I am forcing the same problem. I would also like to implement a non-Newtonian transport model for the solver twoPhaseEulerFoam.



However, I am at the beginning of OpenFOAM programming an I want to generally ask what are the major implementation steps to solve this problem? Which files and folders have to be compiled?


Thanks in advance and best regards,
hebe