[hcl734]

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

[olivierG]

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

[hcl734]

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

[hcl734]

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

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

And an entry in createFluidsField.H

    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

        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

{
    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.

[hcl734]

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)

[olivierG]

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

[hcl734]

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.

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       polynomial;
    thermo          hPolynomial;
    equationOfState icoPolynomial;
    specie          specie;
    energy          sensibleEnthalpy;
}

with chtMultiRegionSimpleFoam



When I look at heRhoThermo.C

#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!

[olivierG]

hello,

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

regards,
olivier

[hcl734]

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.

[hcl734]

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