[jdw135]

I'm working on a problem that requires me to modify the heat source in the flow field. For instance, I am wondering how to change the rho*r term in Eqn. 15 at this page:

http://https://cfd.direct/openfoam/energy-equation/

It says the term, rho*r, is a heat source, where r is any specific heat source. I'd like to make this heat source a function of pressure everywhere in the domain. Is this possible with a fvOptions code? I've seen fvOptions examples using a temperature hot spot, but never one that includes an actual term for the energy equation.

I would like to do this in rhoPimpleFoam if possible, where the energy equation is:

(	fvm::ddt(rho, he)+fvm::div(phi, he)+fvc::ddt(rho, K)+fvc::div(phi, K)+(he.name()=="e" ? fvc::div(fvc::absolute(phi/fvc::interpolate(rho), U), p, "div(phiv,p)") :-dpdt) - fvm::laplacian(turbulence->alphaEff(), he) 	= =fvOptions(rho, he)	)

I'm assuming the fvOptions portion at the end is the part that would allow me to do such a thing? Is anyone knowledgeable about this? I'm fairly new to OpenFOAM and this seems incredibly complicated.

[HPE]

Yes, you can create your own fvOptions to affect this governing equation through the "fvOptions(rho, he)" term.

Please have a look at the available fvOptions here: https://develop.openfoam.com/Develop...ources/derived .

Find the most similar one to your case, and try to modify it accordingly.

[jdw135]

Originally Posted by HPE

Yes, you can create your own fvOptions to affect this governing equation through the "fvOptions(rho, he)" term.

Please have a look at the available fvOptions here: https://develop.openfoam.com/Develop...ources/derived .

Find the most similar one to your case, and try to modify it accordingly.

Thanks for the link! That is what I was looking for. I think the buoyancyEnergy is closest to what I need, as it applies an equation to the entire domain, however, I'm not sure how to include other fields.

For instance, my equation is an isentropic relationship between pressure and heat, so how would I include the mach number and temperature? What about using the pressure from the current and previous time steps? Does anyone have experience with this?

[HPE]

Might help? There are a bit more involved fvOptions here, where temperature field has been used as well: https://develop.openfoam.com/Develop...dels/fvOptions

[jdw135]

Originally Posted by HPE

Might help? There are a bit more involved fvOptions here, where temperature field has been used as well: https://develop.openfoam.com/Develop...dels/fvOptions

Those actually help quite a bit. I'm essentially using as many examples as I can to fill in the parts that I need. Having said that, I need to apply a minimum value to the field I compute. Is there some ability to do something like

if (field < 0.2)
{
field = 0.2;
}

to replace only values within the newly-computed field that meet the if statement? I know I can use a nested for loop, but I doubt that's the correct way to do it.

[HPE]

Hi,
- maybe `Foam::volScalarField& Foam::bound(volScalarField& vsf, const dimensionedScalar& lowerBound)`, like in `bound(k_, this->kMin_);` of `kEpsilon.C`?
- not the second parameter is not `scalar` but `dimensionedScalar`.

[jdw135]

Should it still be a dimensionedScalar if I'm computing a non-dimensional quantity? My thought is to do mach number like this. From looking at other files, it seems like this goes in the constructor in the .c file, so I was planning to implement it like this:

option(sourceName, modelType, dict, mesh),
    UName_(coeffs_.getOrDefault<word>("U", "U")),
    a0_
    (
        dimensionedScalar
        (
            dimLength/dimTime,
            coeffs_.getCheckOrDefault<scalar>
            (
                "a0",
                329.204,
                scalarMinMax::ge(SMALL)    // do I need this line?
            )
        )
    ),
    mach_
    (
        dimensionedScalar
        (
          coeffs_.getCheckOrDefault<scalar>
          (
                "mach",
                U/a0,    // can I put an eqn in the constructor?
                <bound mach number here?>
          )
        )
    )

I feel like I probably can't do the equation in the constructor, as that seems to be setting the default value. Would I create a member function to compute the mach number, or would that go in the system/fvOptions file?

[HPE]

Should it still be a dimensionedScalar if I'm computing a non-dimensional quantity?

Yes. e.g. `dimensionedScalar(dimless, Zero);` or `dimensionedScalar("someName", dimless, 1.0);`

I don't think you can compile the `mach_` like the above.

Could you please explain us (but like explaining to a person with a two-digit IQ) what you want to do as a pseudo-code, if possible?

Many thanks

[jdw135]

Originally Posted by HPE

Yes. e.g. `dimensionedScalar(dimless, Zero);` or `dimensionedScalar("someName", dimless, 1.0);`

I don't think you can compile the `mach_` like the above.

Could you please explain us (but like explaining to a person with a two-digit IQ) what you want to do as a pseudo-code, if possible?

Many thanks

My goal is to add a source term to the rhoPimpleFoam energy equation using fvOptions. This source term is an isentropic relation between pressure and heat, in the form



I was thinking that I would add this by doing the following:

1. compute speed of sound in the flow field
2. compute mach number in the flow field
3. compute the above equation

I'm unsure how to do this, because I have EXTREMELY little experience with OpenFOAM (I've run some tutorials and I have a rhoPimpleFoam case of my flow field that runs). After some reading online, I am sure that using fvOptions is the correct way forward, but I don't have the experience to modify existing options without asking a lot of questions. I'm using a lot of the existing fvOptions codes to help guide me, but the process is slow and painful. I'm still unsure if my approach is reasonable.

[mAlletto]

Hello,


just out of interest, can you provide a reference for the equation.

[jdw135]

Originally Posted by mAlletto

Hello,


just out of interest, can you provide a reference for the equation.

Shapiro, A. (1953) The Dynamics and Thermodynamics of Compressible Fluid Flow. McGraw-Hill, New York.

Specifically, this is a modified derivation from Table 8.1 in that book.

[HPE]

Let's see what we can do during the next week. (At least something for you to progress further rather than a complete solution.)

[mAlletto]

You can use a coded fvOption to include your source term. https://www.openfoam.com/documentati...ces-coded.html

And you also need to decide your equation by dt to insert it into the energy equation. Did you check the dimensions of the source term in the energy equation

[jdw135]

Originally Posted by mAlletto

You can use a coded fvOption to include your source term. https://www.openfoam.com/documentati...ces-coded.html

And you also need to decide your equation by dt to insert it into the energy equation. Did you check the dimensions of the source term in the energy equation

I've been looking at the units, yes. I just posted a generic version for the sake of discussion. I'm working on coding a custom fvOptions right now, but my experience in C++ and OpenFOAM is low, so it's quite a challenge so far. Also, there are a lot of fvOptions files, but many of them do the same task in different ways, so it's hard to figure out what commands are doing what. I have put together a starting point for what I think may work, if you or anyone else is interested in helping me debug it. I would be forever grateful for any help.

[mAlletto]

Maybe it's better if you post the final equation you want to implement. This would lead you faster to your goal.

[jdw135]

Originally Posted by mAlletto

Maybe it's better if you post the final equation you want to implement. This would lead you faster to your goal.

From a quick check of a single term in the energy equation,



I see that the units are



and



So, the units for the first term, and therefore all terms, is



correct? My equation, modified using perfect gas assumption, is



for which the units are



So, if I were to insert my equation into the openfoam energy equation as



it looks like the units are correct, unless I've made a mistake somewhere.

[jdw135]

Here are my first attempts at modifying the variableHeatTransfer fvOptions files to do what I am trying to do. I'm sure there are mistakes, and probably a lot of places where I have misunderstood what the code is doing. To anyone willing to look it over and check it out, I'd be extremely thankful.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2015 OpenFOAM Foundation
    Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::fv::variableHeatTransfer

Group
    grpFvOptionsSources

Description
    Applies an additional source term to the rhoPimpleFoam energy equation via
    the fvOptions.  The source term is computed using Shapiro's influence
    coefficient table as an isentropic relationship between heat and Static
    Pressure, where

    \vartable
      a0      | Speed of Sound
      Cp0     | Specific Heat Capacity (constant pressure)
      mach    | Mach Number
    \endvartable

------------------- this section requires modification -------------------------
    Required fields:
    \verbatim

    \endverbatim

Usage
    Minimal example by using \c constant/fvOptions:
    \verbatim
    variableHeatTransfer1
    {
        // Mandatory entries (unmodifiable)
        type             variableHeatTransfer;

        // Optional entries (runtime modifiable)
        UNbr             U;
        PNbr             p;
        rhoNbr           rho;
        Cp               1005;
        R                287.058;
        a0               329.204;

        // Mandatory/Optional (inherited) entries
        ...
    }
    \endverbatim

    where the entries mean:
    \table
      Property  | Description                     | Type   | Reqd | Dflt
      type      | Type name: variableHeatTransfer | word   | yes  | -
      UNbr      | Name of velocity field          | word   | no   | U
      PNbr      | Name of pressure field          | word   | no   | p
      rhoNbr    | Name of density field           | word   | no   | rho
      Cp        | Specific Heat Capacity          | scalar | no   | 0
      R         | Universal Gas Constant          | scalar | no   | 0
      a0        | Speed of Sound                  | scalar | no   | 0
    \endtable

See also
  - Foam::fv::interRegionHeatTransferModel
  - Foam::fv::tabulatedHeatTransfer
  - Foam::fv::tabulatedNTUHeatTransfer
  - Foam::fv::constantHeatTransfer

SourceFiles
    variableHeatTransfer.C

\*---------------------------------------------------------------------------*/

#ifndef variableHeatTransfer_H
#define variableHeatTransfer_H

#include "fvOption.H"
#include "autoPtr.H"

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

namespace Foam
{
namespace fv
{

/*---------------------------------------------------------------------------*\
                    Class variableHeatTransfer Declaration
\*---------------------------------------------------------------------------*/

class variableHeatTransfer
:
    public option
{
    // Private Data

        //- Name of operand neighbour velocity field
        word UName_;

        //- Speed of Sound
        dimensionedScalar a0_;

        //- Specific Heat Capacity
        const dimensionedScalar Cp0_;

        //- Mach Number
        autoPtr<volScalarField> mach_;

        //- Universal Gas Constant
        const dimensionedScalar R_;


public:

    //- Runtime type information
    TypeName("variableHeatTransfer");


    // Constructors

        //- Construct from dictionary
        variableHeatTransfer
        (
            const word& name,
            const word& modelType,
            const dictionary& dict,
            const fvMesh& mesh
        );

        //- No copy construct
        variableHeatTransfer(const variableHeatTransfer&) = delete;

        //- No copy assignment
        void operator=(const variableHeatTransfer&) = delete;


    //- Destructor
    virtual ~variableHeatTransfer() = default;

    // Member Functions

        //- Compute the Mach Number for the entire field
        virtual void heatCompute
        (
          const dimensionedScalar& a0,
          const dimensionedScalar& Cp0,
          const dimensionedScalar timeStepSize_,
          const volScalarField& rho,
          fvMatrix<scalar>& eqn,
          const label fieldi
        );

    // Public Functions

        //- Read dictionary
        virtual bool read(const dictionary& dict);
};


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

} // End namespace fv
} // End namespace Foam

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

#endif

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

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2016 OpenFOAM Foundation
    Copyright (C) 2020 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "variableHeatTransfer.H"
#include "turbulentFluidThermoModel.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
namespace fv
{
    defineTypeNameAndDebug(variableHeatTransfer, 0);
    addToRunTimeSelectionTable(option, variableHeatTransfer, dictionary);
}
}


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

Foam::fv::variableHeatTransfer::variableHeatTransfer
(
    const word& name,
    const word& modelType,
    const dictionary& dict,
    const fvMesh& mesh
)
:
    option(sourceName, modelType, dict, mesh),
    UName_(coeffs_.getOrDefault<word>("UNbr", "U")),
    PName_(coeffs_.getOrDefault<word>("PNbr","p")),
    rhoName_(coeffs_.getOrDefault<word>("rhoNbr","rho")),
    Cp0_
    (
        dimensionedScalar
        (
            sqr(dimLength)/sqr(dimTime)/dimTemperature,
            coeffs_.getCheckOrDefault<scalar>
            (
                "Cp0",
                1005.0,
                scalarMinMax::ge(SMALL)
            )
        )
    ),
    R_
    (
        dimensionedScalar
        (
            sqr(dimLength)/sqr(dimTime)/dimTemperature,
            coeffs_.getCheckOrDefault<scalar>
            (
                "R",
                287.058,
                scalarMinMax::ge(SMALL)
            )
        )
    ),
    a0_
    (
        dimensionedScalar
        (
            dimLength/dimTime,
            coeffs_.getCheckOrDefault<scalar>
            (
                "a0",
                329.204,
                scalarMinMax::ge(SMALL)
            )
        )
    ),
{

}


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

void Foam::fv::variableHeatTransfer::heatCompute()
{
    const auto& nbrMesh = mesh_.time().lookupObject<fvMesh>(nbrRegionName());

    const auto& UNbr = nbrMesh.lookupObject<volVectorField>(UName_);

    const auto& PNbr = nbrMesh.lookupObject<volScalarField>(PName_);

    const auto& rhoNbr = nbrMesh.lookupObject<volScalarField>(rhoName_);

    tmp<volScalarField> Mach(mag(UNbr)/a0_);

    // Need to bound the lower limit of mach number here...

    eqn += fvm::ddt(((PNbr*timeStepSize_)*(1-(Mach*Mach))*Cp0_)/(1.4*(Mach*Mach)*rho*R_));

}


bool Foam::fv::variableHeatTransfer::read(const dictionary& dict)
{
    if (interRegionHeatTransferModel::read(dict))
    {
        coeffs_.readIfPresent("UNbr", UName_);

        coeffs_.readIfPresent("a0", a0_);
        coeffs_.readIfPresent("Cp0", Cp0_);
        coeffs_.readIfPresent("M", mach_);
        coeffs_.readIfPresent("R", R_);

        return true;
    }

    return false;
}


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

[mAlletto]

Did you try to compile your fvOption and execute it within your solver?

[jdw135]

I haven't yet, mostly because I was trying to find out how to properly bound the mach number before starting on that.

I also am not sure how to compile these, since I've never done it. Do I have to rebuild all of OpenFOAM or just the files I changed?

[mAlletto]

Yes you need to compile only your fvOption. You'll find plenty of tutorials by googling how to write something new by Coping an existing file in OF