CFD Online Discussion Forums - thermophysicalProperties with temperature dependance

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - thermophysicalProperties with temperature dependance (https://www.cfd-online.com/Forums/openfoam-solving/164159-thermophysicalproperties-temperature-dependance.html)

thermophysicalProperties with temperature dependance

Hi all,

is it possible to set upper and lower limits for let's say density when using a polynomial approach?
I want to establish something like that:

Density as a function of temperature
rho_min: 1000
T < 240: 1100
T>240 : C1* T^2 + C2 *T + C3

It would also be possible to use step functions, if available.

Best

hello,
This is not a feature in OF, but it's easy to add this : check how to modify the polynomial density, and add your limiter.
Take a loot at :http://www.tfd.chalmers.se/~hani/kur...nFoam%20v2.pdf

This is for viscosity, but do the same for density...
The basic step are:
1) copy $FOAM_SRC/thermophysicalModels/specie/equationOfState/icoPolynomial/ where you want.
2) rename all "icoPolynomial" by what you want like "icoLimitPolynomail"(file name AND inside !)
3) modify in file to add what you want
4) add a make file (make sur lib will be in FOAM_USER_LIBBIN)
4) compile (wmake) and use it ! (do not forget the libname in you controlDict)

regards,
olivier

Ok found the solution

Use icoPolynomial for density as explained here http://cfd.direct/openfoam/user-guide/thermophysical/

and set the upper and lower limits for rho in the fvSolution

Thanks for showing me that paper, but to be honest I'm not very familiar with OF programming and I found this simple solution for density which works good for me.

BUT I want to do the same thing for heat capacit in chtMultiRegionSimpleFoam-solver.
For density the feature is implemented in the pressure equation as

Quote:

rho = thermo.rho();
rho = max(rho, rhoMin[i]);
rho = min(rho, rhoMax[i]);
rho.relax();

And an entry in createFluidsField.H

Code:



    PtrList<dimensionedScalar> CpMax(fluidRegions.size());

    PtrList<dimensionedScalar> CpMin(fluidRegions.size());





        rhoMax.set(i, new dimensionedScalar(simpleDict.lookup("rhoMax")));

        rhoMin.set(i, new dimensionedScalar(simpleDict.lookup("rhoMin")));

And added in createFluidsField.H

Code:

        CpMax.set(i, new dimensionedScalar(simpleDict.lookup("CpMax")));

        CpMin.set(i, new dimensionedScalar(simpleDict.lookup("CpMin")));

But when I take a look at the EEqn.H

Code:

{

    volScalarField& he = thermo.he();



    fvScalarMatrix EEqn

    (

        fvm::div(phi, he)

      + (

            he.name() == "e"

          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))

          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))

        )

      - fvm::laplacian(turb.alphaEff(), he)

     ==

        rad.Sh(thermo)

      + fvOptions(rho, he)

    );



    EEqn.relax();



    fvOptions.constrain(EEqn);



    EEqn.solve();



    fvOptions.correct(he);



    thermo.correct();

    rad.correct();



    Info<< "Min/max T:" << min(thermo.T()).value() << ' '

        << max(thermo.T()).value() << endl;

}

I can't even see where they are using specific heat capacity?
Could you help me with this one?
It would also be nice to write the Cp-field while solving to check whether it's correct or not.

Viscosity depending on temperature exponentially

2) Question

Is it also possible to set an exponential dependance between viscosity and temperature?
Something like:

(C1 + C2*exp(-(T-Tref)/C5) + C3*exp(-(T-Tref)/C6) + C4*exp(-(T-Tref)/C7)

hello,

You are wrong about the limit in fvSolution => this bound rho, but NOT in a conservative way. This is only for stability purpose. You should use a correct model for density.

About 1°): Cp is inside alphaEff() ( ~ something like k/Cp + kt/Cp)
You can use polynomal for Cp => use janaf model

About 2°) by defaut, OF has only Sutherland viscosity(T) for gas. Take a look at the link i give you: this is how to add a custom viscosity.

BTW, you should not be afraid about coding with OF: this is the way do do.

regards,
olivier

Thanks for your advice, I followed the tutorial you showed me for implementing a temperature dependent viscosity, but this tutorial is only for incompressible flow using a viscosityProperties dictionary. But I want to calculate compressible using thermophysicalProperties.

Code:

thermoType

{

    type            heRhoThermo;

    mixture         pureMixture;

    transport       polynomial;

    thermo          hPolynomial;

    equationOfState icoPolynomial;

    specie          specie;

    energy          sensibleEnthalpy;

}

with chtMultiRegionSimpleFoam

When I look at heRhoThermo.C

Code:



#include "heRhoThermo.H"



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



template<class BasicPsiThermo, class MixtureType>

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

{

    const scalarField& hCells = this->he().internalField();

    const scalarField& pCells = this->p_.internalField();



    scalarField& TCells = this->T_.internalField();

    scalarField& psiCells = this->psi_.internalField();

    scalarField& rhoCells = this->rho_.internalField();

    scalarField& muCells = this->mu_.internalField();

    scalarField& alphaCells = this->alpha_.internalField();



    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]);

        alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);

    }



    forAll(this->T_.boundaryField(), patchi)

    {

        fvPatchScalarField& pp = this->p_.boundaryField()[patchi];

        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];

        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];

        fvPatchScalarField& prho = this->rho_.boundaryField()[patchi];



        fvPatchScalarField& ph = this->he().boundaryField()[patchi];



        fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];

        fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];



        if (pT.fixesValue())

        {

            forAll(pT, facei)

            {

                const typename MixtureType::thermoType& mixture_ =

                    this->patchFaceMixture(patchi, facei);



                ph[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]);

                palpha[facei] = mixture_.alphah(pp[facei], pT[facei]);

            }

        }

        else

        {

            forAll(pT, facei)

            {

                const typename MixtureType::thermoType& mixture_ =

                    this->patchFaceMixture(patchi, facei);



                pT[facei] = mixture_.THE(ph[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]);

                palpha[facei] = mixture_.alphah(pp[facei], pT[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)

    {

        Info<< "entering heRhoThermo<MixtureType>::correct()" << endl;

    }



    calculate();



    if (debug)

    {

        Info<< "exiting heRhoThermo<MixtureType>::correct()" << endl;

    }

}





// ************************************************************************* //

I can't see much resemblance between this and the viscosityModel files.
Where should I start when modifying the compressible case?
I mean basically I only need an exponential function for defining mu and I need to know how to write the transportProperties as fields.
Thx again for your help!

hello,

You should add a custom thermo model.
Check this : http://www.cfd-online.com/Forums/ope...tml#post451937

regards,
olivier

Ok I checked the post you recommended to me and I am trying to get into programming with OpenFoam.
But since this is quite tedious I am wondering whether it is possible to define a kind of lookup table for a thermophysical propertie in OpenFoam where I can put some data points (calculated externally) and OpenFoam does linear interpolation between those.
I hope you get my idea.

Didn't checked it but for everybody coming after me

This gives you the possibility to assign exponential terms
https://github.com/OpenFOAM/OpenFOAM...1ec9af6f399210