[EnricoA]

Hello alexeym,

I am trying to implement a source term for turbulent dissipation in a simulation using the standard k-epsilon model following your explanation about the scalarCodedSource. I ended up with this fvOptions file (where there is also a source for the momentum equation):

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


momentumSource
{
    type            vectorSemiImplicitSource;
    selectionMode   cellZone;
    cellZone        disk;
    active          true;

    vectorSemiImplicitSourceCoeffs
    {
        volumeMode      absolute;
        injectionRateSuSp
        {
            U           ((0 0 -180e3) 0);
        }
    }
}

dissipationSource
{
    type            scalarCodedSource;
    selectionMode   cellZone;
    cellZone        disk;
    active          true;

    scalarCodedSourceCoeffs
    {
        fieldNames      (epsilon);
        redirectType    sourceTime;

        codeInclude
        #{

        #};

        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};

        codeAddSup
        #{
            const Time& time = mesh().time();
            const scalarField& V = mesh_.V();
            const vectorField& C = mesh_.C();
            const scalarField& nut_ = nut();
            const vectorField& U_ = U();
            const scalarField& k_ = k();
            const scalarField& Pt = nut_*2*magSqr(symm(fvc::grad(U_)));
            scalarField& epsilonSource = eqn.source();
            forAll(C, i)
            {
                epsilonSource[i] -= 0.37*Pt[i]*Pt[i]/k_[i]*V[i];
            }
            Pout << "***codeAddSup***" << endl;
        #};

        codeSetValue
        #{
                Pout<< "**codeSetValue**" << endl;
        #};

        // Dummy entry. Make dependent on above to trigger recompilation
        code
        #{
            $codeInclude
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
    }

    sourceTimeCoeffs
    {
        // Dummy entry
    }
}

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

Apparently there are some mistakes in my file which I cannot find, because when I run the simulation the output is:

--> FOAM Warning : 
From function void option::checkApplied() const
in file fvOptions/fvOption.C at line 368
Source dissipationSource defined for field epsilon but never used

Do you have any advises to give me?

Thanks

[alexeym]

Hi,

My advice: implement your own turbulence model I.e. you take k-epsilon family model, copy it, rename it, add your code.

There are certain conditions for fvOptions to work, see, for example, UEqn.H of pimpleFoam:

tmp<fvVectorMatrix> UEqn
(
    fvm::ddt(U)
  + fvm::div(phi, U)
  + turbulence->divDevReff(U)
 ==
    fvOptions(U)
);

UEqn().relax();

fvOptions.constrain(UEqn());

See all these calls to fvOptions. Now if you look at kEpsilon.C, there are no such calls.

During last workshop there was an idea of implementation of turbulence models, which use fvOptions framework. I do not know if there is any progress, for me it is still on TODO list.

[sandeeprapol]

Hello everyone
I am implementing volumetric heat source in chtmultiregionsimplefoam with help of fvoptions files
in that "duration ....." and " h(.... 0)" so my problem is

1) how to convert 200 W source into enthalpy h(... 0) with help of "duration ...." in sec
2) in that "fvOptions" file "duration..." so which time put here in duration

[zfaraday]

Originally Posted by sandeeprapol

Hello everyone
I am implementing volumetric heat source in chtmultiregionsimplefoam with help of fvoptions files
in that "duration ....." and " h(.... 0)" so my problem is

1) how to convert 200 W source into enthalpy h(... 0) with help of "duration ...." in sec
2) in that "fvOptions" file "duration..." so which time put here in duration

Hello Sandeep!

1) You don't need to convert nothing into nothing. For the case of using chtMultiRegionSimpleFoam you just have to put the thermal power value into the first place within the brackets, like that:

energySource
{
    type            scalarSemiImplicitSource;
    active          true;
    selectionMode   all;

    scalarSemiImplicitSourceCoeffs
    {
        volumeMode      absolute;//specific;//
        injectionRateSuSp
        {
            h           (q 0); //   q in [W]; or in [W/m³] if you use specific mode
        }
    }
}

2) Never used the duration field (I don't even know wether it exist) as you can see in the piece of code above.

Hope it helps.

Best regards,

Alex

[alexeym]

@faraday

Originally Posted by zfaraday

2) Never used the duration field (I don't even know wether it exist) as you can see in the piece of code above.

In fact the field comes from base abstract fvOption and defines duration of time the option is active:

inline bool Foam::fv::option::inTimeLimits(const scalar time) const
{
    return
    (
        (timeStart_ < 0)
     ||
        (
            (mesh_.time().value() >= timeStart_)
         && (mesh_.time().value() <= (timeStart_ + duration_))
        )
    );
}

bool Foam::fv::option::isActive()
{
    if (active_ && inTimeLimits(mesh_.time().value()))
    {
        ---
        return true;
    }
    else
    {
        return false;
    }
}

[zfaraday]

Thanks for the clarification Alexey! I didn't remember that because I never had to specify the duration of the source.

Everything must be clear now for @Sandeep.

[sandeeprapol]

hello alex,
thank you for the replay , I used that syntax in fvOptions file but result showing no heat generation, it shows constant temperature

when I use my fvOptins that showing heat generation it include "duration....." of time
-------#---------------------------#-----------------------------#--------------------
heatSource
{
type scalarSemiImplicitSource;
active on;
timeStart 0.;
duration 1e3;
selectionMode cellSet;
cellSet IC1;

scalarSemiImplicitSourceCoeffs
{
// volumeMode absolute; // Values are given as <quantity>
volumeMode specific; // Values are given as <quantity>/m3

injectionRateSuSp // Semi-implicit source term S(x) = S_u + S_p x
{
h (200000 0);
}
}
}
--------------------#--------------------------------#-----------------------#------------------#---
am trying to implement 200 W/m3 volumetric heat generating source for circuit board cooling

[zfaraday]

Well, I don't understand why my specification is not working. As per what I see you are using a cellSet as a selection mode, while I sellect all cells belonging to one region. If you just copypasted all my specification of course it's not going to work...

Another point I don't get is why you are setting a generation of 200000 when you say your generation is supposed to be of 200 W/m2... If you give a higher value than the one you need, obviously the effect will be more visual and bigger...

Regards,

Alex

Ps: excuse me if I wrote something wrong, I'm writing from my phone

[sandeeprapol]

hello alex,
you are on the right way you see in src/fvOptios/lninclude/fvOptionListTemplates.C
at line no 135
ds = rho.dimensions()*fld.dimensions()/dimTime*dimVolume
rho * h *1/time * vol
kg/m3 * j/kg * m3 / sec
j/s
W
so here is dimTime in sec plese tell me how to implement 200 W/m3 heat generating source

[zfaraday]

check this out!

[sandeeprapol]

hello alex
thank you for your humble reolay
in fvoptions
volumeMode absolute; // Values are given as <quantity>
volumeMode specific; // Values are given as <quantity>/m3
what is "absolute" and "specific"
if my heat generating element is 200 W having dimension 0.2m*0.01m*0.05m
1) if I consider "absolute" then total element given value is 200 W is this correct or not?
2) if I consider "specific" then total element each cell point given value is 200 W is this correct or not?

[zfaraday]

Dear Sandeep,

the meaning of "absolute" and "specific" is given behind "//". Then, if you consider "absolute" mode you just have to give the constant source term a value of 200. On the other hand, if you prefer to use "specific" mode you have to use a value of 200/(0.2*0.01*0.05).

That's all!

Best regards,

Alex

[sandeeprapol]

thank you alex,
which parameter analyze in post processing.I thought that only temp of oullet, & velocity is analyze in chtmultiregionSimpleFoam and chtmultiRegionFoam please suggest me

[EnricoA]

Hi,

I would like to follow up my previous post regarding the implementation of a source term for the turbulent dissipation rate with scalarCodedSource.
I eventually implemented my own turbulence model of the k-epsilon family. This is the main modification which takes into account the possible source term:

// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
    fvm::ddt(epsilon_)
  + fvm::div(phi_, epsilon_)
  - fvm::laplacian(DepsilonEff(), epsilon_)
 ==
    C1_*G*epsilon_/k_
  - fvm::Sp(C2_*epsilon_/k_, epsilon_)
  + fvOptions(epsilon_)
);

epsEqn().relax();
fvOptions.constrain(epsEqn());
epsEqn().boundaryManipulate(epsilon_.boundaryField());
solve(epsEqn);
fvOptions.correct(epsilon_);
bound(epsilon_, epsilonMin_);

Then I was trying to add a source of the following form in a particular region:

S = 0.37 G^2 / k

where G is the the generation of turbulent kinetic energy computed in the sandard k-epsilon model and k is the turbulent kinetic energy. In order to do this, I added the following piece of code to fvOptions:

dissipationSource
{
    type            scalarCodedSource;
    selectionMode   cellZone;
    cellZone        diss;
    active          true;

    scalarCodedSourceCoeffs
    {
        fieldNames      (epsilon);
        redirectType    sourceTime;

        codeInclude
        #{
            #include "fvCFD.H"
        #};

        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};

        codeAddSup
        #{
            const scalarField& V = mesh_.V();
            const vectorField& C = mesh_.C();
            const volScalarField& nut = mesh().lookupObject<volScalarField>("nut");
            const volVectorField& U = mesh().lookupObject<volVectorField>("U");
            const volScalarField& k = mesh().lookupObject<volScalarField>("k");
            const volScalarField& G = nut*2*magSqr(symm(fvc::grad(U))); // as computed in the k-epsilon model
            scalarField& epsilonSource = eqn.source();
            forAll(C,i)
            {
                epsilonSource[i] -= 0.37*G[i]*G[i]/k[i];
            }
            Pout << "***codeAddSup***" << endl;
        #};

        codeSetValue
        #{
            Pout<< "**codeSetValue**" << endl;
        #};

        // Dummy entry. Make dependent on above to trigger recompilation
        code
        #{
            $codeInclude
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
    }

    sourceTimeCoeffs
    {
        // Dummy entry
    }
}

When I run a simulation, the compilation of the dynamic library doesn't give errors but apparently there is something wrong with my implementation, since the results that I get don't match with the expected one.
Can you see any errors in the code?

Thanks

[nimasam]

Dear Alex

How to set multiple discrete heat sources using scalarSemiImplicitSource ?

[student666]

Originally Posted by zfaraday

{
type scalarSemiImplicitSource;
active true;
selectionMode all;

scalarSemiImplicitSourceCoeffs
{
volumeMode absolute;//specific;//
injectionRateSuSp
{
h (q 0); // q in [W]; or in [W/m³] if you use specific mode
}
}
}

Hi,

just a clarification for buoyantSimpleFoam.
This solver uses EEqn.C file, so you have to specify the thermal model.
In case you set the following dict for thermophysicalProperties

    {
        type heRhoThermo;
        mixture pureMixture;
        transport const;
        thermo hConst;
        equationOfState perfectGas;
        specie specie;
        energy sensibleEnthalpy;
    }

q parameter (specific) is equal to Watt only if you consider a fluid with Cp = 1kJ/(kg*K) and rho = 1 kg/m^3 as for Q= rho cp * T(1Kelvin)/1sec = 1W/m3 as h = Cp*T(1Kelvin)/1sec = 1W/kg

but if you consider a fluid with Cp = 4.186kJ/(kg*K) and rho = 1000 kg/m3, if you want to have 1W/m3 you have to set h value (q parameter [specific] ) equal to 0,000000239

 q = Q = rho*Cp*T/t ==> 1/1000/4186

is it correct?

[student666]

I performed this simple test case.
3D simulation.
All patches have been set to "walls"(T=293K), but patch named "fondo" has been set to zeroGradient for T.
cellZone "c0" , at the bottom (coord z=0), has 500W power source.

laminar.

simulation performed with openFoam v4.1

run

blockMesh
topoSet
renumberMesh -overwrite
buoyantPimpleFoam > log.1

when steadyness is reached, stop running and

wallHeatFlux

you may see dissipation through patch "pareti" equal to heat (absolute) source 500W in fvOptions.

Time = 2.554868161131
Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

Reading/calculating face flux field phi

Selecting turbulence model type laminar

Wall heat fluxes [W]
pareti
    convective: -497.79226
    radiative:  -0
    total:      -497.79226
fondo
    convective: 0
    radiative:  -0
    total:      0

Time = 2.9167083871732
Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

Reading/calculating face flux field phi

Selecting turbulence model type laminar

Wall heat fluxes [W]
pareti
    convective: -498.87276
    radiative:  -0
    total:      -498.87276
fondo
    convective: 0
    radiative:  -0
    total:      0

Time = 3.27661601896623
Selecting thermodynamics package 
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

Reading/calculating face flux field phi

Selecting turbulence model type laminar

Wall heat fluxes [W]
pareti
    convective: -499.42118
    radiative:  -0
    total:      -499.42118
fondo
    convective: 0
    radiative:  -0
    total:      0

End

Hope this can clarify any further doubt...mine too.

Regards.

[ChrisBa]

Hello,

I'm trying to implement a volumetric source with scalarCodedSource.
At the surface the equation: q=q0 * e^(re^2/r^2)
In the depth the value q0 is described with functions.

I'm using Openfoam 4.1 and the solver buoyantBoussinesqSimpleFoam.

After calling the function I get this:

Creating finite volume options from "constant/fvOptions"



--> FOAM FATAL ERROR: 
Attempt to return primitive entry ITstream : /home/cfd/OpenFOAM/cfd-4.1/run/elektronen/constant/fvOptions.heatSource.scalarCodedSourceCoeffs.codeInclude, line 31, IOstream: Version 2.0, format ASCII, line 0, OPENED, GOOD
    primitiveEntry 'codeInclude' comprises 
        on line 31 the verbatim string "\
           \
        "
 as a sub-dictionary

    From function virtual const Foam::dictionary& Foam::primitiveEntry::dict() const
    in file db/dictionary/primitiveEntry/primitiveEntry.C at line 189.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::primitiveEntry::dict() const at primitiveEntry.C:?
#3  Foam::dictionary::substituteScopedKeyword(Foam::word const&) at ??:?
#4  Foam::entry::New(Foam::dictionary&, Foam::Istream&) at ??:?
#5  Foam::dictionary::read(Foam::Istream&, bool) at ??:?
#6  Foam::dictionary::dictionary(Foam::fileName const&, Foam::dictionary const&, Foam::Istream&) at ??:?
#7  Foam::dictionaryEntry::dictionaryEntry(Foam::keyType const&, Foam::dictionary const&, Foam::Istream&) at ??:?
#8  Foam::entry::New(Foam::dictionary&, Foam::Istream&) at ??:?
#9  Foam::dictionary::read(Foam::Istream&, bool) at ??:?
#10  Foam::dictionary::dictionary(Foam::fileName const&, Foam::dictionary const&, Foam::Istream&) at ??:?
#11  Foam::dictionaryEntry::dictionaryEntry(Foam::keyType const&, Foam::dictionary const&, Foam::Istream&) at ??:?
#12  Foam::entry::New(Foam::dictionary&, Foam::Istream&) at ??:?
#13  Foam::dictionary::read(Foam::Istream&, bool) at ??:?
#14  Foam::dictionary::dictionary(Foam::fileName const&, Foam::dictionary const&, Foam::Istream&) at ??:?
#15  Foam::dictionaryEntry::dictionaryEntry(Foam::keyType const&, Foam::dictionary const&, Foam::Istream&) at ??:?
#16  Foam::entry::New(Foam::dictionary&, Foam::Istream&) at ??:?
#17  Foam::dictionary::read(Foam::Istream&, bool) at ??:?
#18  Foam::operator>>(Foam::Istream&, Foam::dictionary&) at ??:?
#19  Foam::IOdictionary::readFile(bool) at ??:?
#20  Foam::IOdictionary::IOdictionary(Foam::IOobject const&) at ??:?
#21  Foam::fv::options::options(Foam::fvMesh const&) at ??:?
#22  Foam::fv::options::New(Foam::fvMesh const&) at ??:?
#23  ? at ??:?
#24  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#25  ? at ??:?
Abgebrochen (Speicherabzug geschrieben)

My fvOption-file looks like this:

heatSource
{
    type            scalarCodedSource;

    active          true;
    selectionMode   all;

    scalarCodedSourceCoeffs
    {
        fieldNames      (T);
        name            sourceTime;

        codeInclude
        #{
           
        #};

        codeCorrect
        #{
            Pout<< "**codeCorrect**" << endl;
        #};

        codeAddSup
        #{
            scalarField& TSource = eqn.source();


            //Values for checking the equations
            const scalar power = 1000;
            const scalar Radius = 0.5;
            const scalar rho = 1000;
            const scalar Cp = 4.19;
            const scalar xCenter = 3;
            const scalar yCenter = 0.01;
            const scalar zCenter = 2;             
            
            // Equations 

            // -Face centers
            const List<point>& cf = p.Cf();
            const scalar xCF = cf[c][0];
            const scalar XF = xCenter - xCF;
            const scalar check+ = xCenter + Radius; 
            const scalar check- = xCenter - Radius;

            if ((xCF < check+) && (xCF > check-))
            {
                const scalar r2 = XF * XF; 
                const scalar re2 = Radius * Radius; 
                const scalar factor = r2/re2; 
                const scalar pre = exp(factor); 
                const scalar qfa_ [c] = power * exp;
                const scalar TzuF_ [c] = (1/(rho*Cp))*qfa   
            }
            else 
            {
                const scalar qfa_ [c] = 0;
            };
            
            // -Depth centers or cell centers
            
            
            
            
            
                       
         

            TSource -= Tzuf; 
            Pout << "***codeAddSup***" << endl;

        #};
 
        codeSetValue
        #{
            Pout<< "**codeSetValue**" << endl;
        #};

        // Dummy entry. Make dependent on above to trigger recompilation
        code
        {
           $codeInclude
           $codeCorrect
           $codeAddSup
           $codeSetValue
        };
    }

    sourceTimeCoeffs
    {
        // Dummy entry
    }
}

Maybe somebody can help me.

regards
Chris

[hanness]

Hi Chris,

I'm not too much into the scalarCodedSource fvOptions but just by looking at the error message it tells you that on line 31 of your code something is going wrong. More precisely it is in the section codeInclude (I suppose that is line 31 when looking at the entiere fvOptions file including its header.
This section supposedly must be filled and cannot be left blank. As I'm not familiar with this tool I can't help you any further but maybe it helps.
Regards
Hannes

[ChrisBa]

Thanks for your reply.
I look into this.