Calculate integral length scale with coded functions

Siassei · March 18, 2019, 14:38

Hello,

atm I had no experience with coded functions. My first try is to produce a scalar field by calculating the integral length scale of an element and divided the result with the element volume. I would be use an RANS simulation with an kOmega turbulence model.

$k \text{ := turbulent kinetic energy}$
$\varepsilon \text{ := energy dissipation rate}$

$\text{integral length scale} = k ^{3/2} / \varepsilon$

$SField_i = \frac{\text{integral length scale} }{ \text{Element Volume}^\frac{1}{3}}$

my first lines of system/controlDict looks like

Code:

functions
{
    calcMeshFactorForLES   
    {
        type                coded;
        libs                ("libutilityFunctionObjects.so");
        enabled             yes;
        writeControl        timeStep;
        code
        #{
            // get turbulent kinetic energy
            const volScalarField& k    (???);
            // get energy dissipation rate
            const volScalarField& eps  (???);
            //
            auto integralLengthScale = FOAM::pow(k, 3.0 / 2.0) / eps;
            //
            auto volume (???);
            auto factor = integralLengthScale / FOAM::pow(volume, 1.0 / 3.0);
        #}
    }
}

but I'm very new of coding with OF. The $FOAM_TUTORIALS also have a very small set of examples. It would be a pleasure someone tells my the way / solution. It would be very helpfull for me. The first steps are the difficultities

Thanks.

SvenBent · March 19, 2019, 08:56

Hi Thomas.

I have tried to finish you code. I have not had the time to check if it works, but hopefully it will. Basically you use the mesh().lookUpObject<Type>("VarName") to get different fields from the object registry. If it cant find the name of the variable, un-comment the INFO line, and run the solver, it will show you whats available in the terminal.

Code:

functions
{
    calcMeshFactorForLES   
    {
        type                coded;
        libs                ("libutilityFunctionObjects.so");
        enabled             yes;
        writeControl        timeStep;
        code
        #{
            // get turbulent kinetic energy
            const volScalarField& k  =  mesh().lookupObject<volScalarField>("k");
            // get energy dissipation rate
            const volScalarField& eps  = mesh().lookupObject<volScalarField>("epsilon");
            //
            auto integralLengthScale = FOAM::pow(k, 3.0 / 2.0) / eps;
            //

            //Info<<mesh().thisbd()<<endl; //print all objects to the terminal
           const scalarField& volume  = mesh_.V(); // Cell volume
            auto factor = integralLengthScale / FOAM::pow(volume, 1.0 / 3.0);
        #}
    }
}

hope it helps.

best regards
Mathias

Siassei · March 19, 2019, 14:31

Hi Mathias,

thank you for the quick answer, but a few question remain

First of all, I found the codedFunctionObject, that helps a lot.

I fixed something

Code:

    calcMeshFactorForLES
    {
        libs        ("libutilityFunctionObjects.so");
        type coded;
        name calculateFactorForLESField;
        codeWrite
        #{
            // get turbulent kinetic energy
            const volScalarField& k  =  mesh().lookupObject<volScalarField>("k");
            // get energy dissipation rate
            const volScalarField& eps  = mesh().lookupObject<volScalarField>("epsilon");
            //
            auto integralLengthScale = Foam::pow(k, 3.0 / 2.0) / eps;
            //
            //Info<<mesh().thisDb()<<endl; //print all objects to the terminal
            const scalarField& volume  = mesh().V(); // Cell volume
            auto factor = integralLengthScale / Foam::pow(volume, 1.0 / 3.0);
            //Info << factor << endl;
            writeToField???
        #};
    }

Wow, Foam looks very nice to write mathematics in code. I hope the pow functions works correct on the field.

If I do not use the k-epsilon turbulence model, how I can calculate k and epsilon?
Where get I the values for the formula variables?
How I add the result to an declared scalar field?

How and where I declare an scalar field? Directory 0?
0/fator1Field

Code:

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      fator1Field;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -3 0 0 0 0];

internalField   uniform 0;


boundaryField
{
    inlet
    {
        type        fixedValue;
        value       $internalField;
    }

    sym1
    {
        type        empty;
    }

    sym2
    {
        type        empty;
    }

    outlet
    {
        type        fixedValue;
        value       $internalField;
    }

    walls
    {
        type        fixedValue;
        value       $internalField;
    }
}

Siassei · March 19, 2019, 15:42

So, the answer of my last question are solved by my self

For someone that has the same problem

Code:

    calcMeshFactorForLES
    {
        libs        ("libutilityFunctionObjects.so");
        type coded;
        name calculateFactorForLESField;
        codeWrite
        #{
            // get turbulent kinetic energy
            const volScalarField& k  =  mesh().lookupObject<volScalarField>("k");
            // get energy dissipation rate
            const volScalarField& eps  = mesh().lookupObject<volScalarField>("epsilon");
            //
            auto integralLengthScale = Foam::pow(k, 3.0 / 2.0) / eps;
            //
            //Info<<mesh().thisDb()<<endl; //print all objects to the terminal
            const scalarField& volume  = mesh().V(); // Cell volume
            auto myFactor = integralLengthScale / Foam::pow(volume, 1.0 / 3.0);

            volScalarField approxCellSizeForLESFactor
            (
                IOobject
                (
                    "approxCellSizeForLESFactor",
                    time_.timeName(),
                    mesh(),
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
                mesh(),
                dimensionedScalar("approxCellSizeForLESFactor", dimless, scalar(0.0))
            );
            approxCellSizeForLESFactor.primitiveFieldRef() = myFactor;
            approxCellSizeForLESFactor().write();
        #};
    }

Last questions are

I put mesh() to IOobject and volScalarField. Is that correct or should I pass something else?
If I do not use the k-epsilon turbulence model, how I can calculate k and epsilon?
Where get I the values for the formula variables?
Results looking correct. Is the use of the Foam:ow function correct? Do the pow function what the first indent is?

SvenBent · March 21, 2019, 03:17

Great work.

Your IOobject and volScalarField approxCellSizeForLESFactor looks correct to me

Regarding the calculation of k and epsilon. Assuming that you are using LES, I think you have to look in the literature and find the link between the variables used in LES and RANS simulations. I am not very familiar with LES theory or simulations, so you might want to ask someone with more experience in this field than me.

The foam:

ow function is basically the same as the pow function used in C++.
Say you have

where

is some scalarField, you would write that as

Code:

foam::pow(x,2.0)

So yes I would say that you are using the function correctly.

Best regards

Mathias

Siassei · March 21, 2019, 13:55

Thanks Matthias.

Quote:

Regarding the calculation of k and epsilon. Assuming that you are using LES, ...

Wrong direction. I'm using RANS and learning to use FOAM.
The calculated function above is a approximated value to asume mesh size for LES calculation. But you can also use the factor something else.

I try to get an universal function that can also handle RANS calculation if k or epsilon is not present. Please correct me, when I wrong. But I think the k & epsilon field in FOAM-RANS is only calculated, when the turbulence model is included the calculation. I think (without test

) calculation with k-omega turbulence model do not include epsilon field.

Is a way to force to calculate k & epsilon unindepend the choosen turbulence model in RANS?
Or, how I get the needed variables to calculate the shit.

Siassei · March 21, 2019, 14:35

I opened a new thread for epsilon and k calculation.

calf.Z · November 22, 2019, 07:47

If I am doing DNS and I want to get k and epsilon after calculation, how should I do?

sumedhjain · December 17, 2020, 08:31

I used your code to find the meshFactor for LES. Unfortunately, the code executes with a minor error but the factor filed is not being written at all.

OF version - v1912
solver - simpleFoam in postProcess mode

Code:

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1912                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : _f3950763fe-20191219 OPENFOAM=1912
Arch   : "LSB;label=32;scalar=64"
Exec   : simpleFoam -postProcess -parallel -latestTime -decomposeParDict system/decomposeParDict.ptscotch
Date   : Dec 17 2020
Time   : 13:19:03
Host   : ip-172-31-5-127
PID    : 40150
I/O    : uncollated
Case   : /fsx/Fsx/export/Fsx/OF64_6_fa350
nProcs : 288
Hosts  :
(
    (ip-172-31-5-127 96)
    (ip-172-31-15-187 96)
    (ip-172-31-8-42 96)
)
Pstream initialized with:
    floatTransfer      : 0
    nProcsSimpleSum    : 0
    commsType          : nonBlocking
    polling iterations : 0
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 1000


SIMPLE: no convergence criteria found. Calculations will run for 500 steps.

--> FOAM Warning :

    From function virtual bool Foam::functionObjects::codedFunctionObject::read(const Foam::dictionary&)
    in file codedFunctionObject/codedFunctionObject.C at line 288
    Reading "/fsx/Fsx/export/Fsx/OF64_6_fa350/system/controlDict.functions.calcMeshFactorForLES"
    No critical "code" prefixed keywords found.
Please check the code documentation for more details.

Using dynamicCode for functionObject calcMeshFactorForLES at line -1 in "/fsx/Fsx/export/Fsx/OF64_6_fa350/system/controlDict.functions.calcMeshFactorForLES"
Time = 1000
Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RAS
Selecting RAS turbulence model kOmegaSST
Selecting patchDistMethod meshWave
RAS
{
  RASModel        kOmegaSST;
    turbulence      on;
    printCoeffs     on;
    alphaK1         0.85;
    alphaK2         1;
    alphaOmega1     0.5;
    alphaOmega2     0.856;
    gamma1          0.555556;
    gamma2          0.44;
    beta1           0.075;
    beta2           0.0828;
    betaStar        0.09;
    a1              0.31;
    b1              1;
    c1              10;
    F3              false;
    decayControl    false;
    kInf            0;
    omegaInf        0;
}

No MRF models present

No finite volume options present

End

Finalising parallel run

ukaya · February 25, 2021, 13:54

It seems like you couldn't copy the function properly to controldict.

Also try to run it without all the keywords, but the necessary ones:
simpleFoam -postProcess -latestTime

March 18, 2019, 14:38	Calculate integral length scale with coded functions	#1
Siassei Member Join Date: Jun 2011 Location: Germany Posts: 59 Rep Power: 14	Hello, atm I had no experience with coded functions. My first try is to produce a scalar field by calculating the integral length scale of an element and divided the result with the element volume. I would be use an RANS simulation with an kOmega turbulence model. $k \text{ := turbulent kinetic energy}$ $\varepsilon \text{ := energy dissipation rate}$ $\text{integral length scale} = k ^{3/2} / \varepsilon$ $SField_i = \frac{\text{integral length scale} }{ \text{Element Volume}^\frac{1}{3}}$ my first lines of system/controlDict looks like Code: functions { calcMeshFactorForLES { type coded; libs ("libutilityFunctionObjects.so"); enabled yes; writeControl timeStep; code #{ // get turbulent kinetic energy const volScalarField& k (???); // get energy dissipation rate const volScalarField& eps (???); // auto integralLengthScale = FOAM::pow(k, 3.0 / 2.0) / eps; // auto volume (???); auto factor = integralLengthScale / FOAM::pow(volume, 1.0 / 3.0); #} } } but I'm very new of coding with OF. The $FOAM_TUTORIALS also have a very small set of examples. It would be a pleasure someone tells my the way / solution. It would be very helpfull for me. The first steps are the difficultities Thanks.

March 21, 2019, 03:17		#5
SvenBent New Member Mathias Poulsen Join Date: Feb 2018 Location: Denmark Posts: 9 Rep Power: 8	Great work. Your IOobject and volScalarField approxCellSizeForLESFactor looks correct to me Regarding the calculation of k and epsilon. Assuming that you are using LES, I think you have to look in the literature and find the link between the variables used in LES and RANS simulations. I am not very familiar with LES theory or simulations, so you might want to ask someone with more experience in this field than me. The foam:ow function is basically the same as the pow function used in C++. Say you have where is some scalarField, you would write that as Code: foam::pow(x,2.0) So yes I would say that you are using the function correctly. Best regards Mathias

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Difference between mixing length and turbulence length scale	MissCFD	CFX	1	May 7, 2015 12:58
Estimation of turbulence length scale	J.Y.Won	OpenFOAM Pre-Processing	0	April 28, 2015 07:16
Integral Length Scale vs Turbulent Length Scale	Scott Nordsen	FLUENT	2	January 20, 2014 01:17
Turbulent length scale setting	realanony87	CFX	5	October 14, 2012 18:03
Turbulent Characteristic Length/Turbulent Length Scale	fluentnoob	FLUENT	2	July 3, 2009 08:40

March 21, 2019, 14:35		#7
Siassei Member Join Date: Jun 2011 Location: Germany Posts: 59 Rep Power: 14	I opened a new thread for epsilon and k calculation.

November 22, 2019, 07:47		#8
calf.Z Senior Member Jianrui Zeng Join Date: May 2018 Location: China Posts: 157 Rep Power: 7	If I am doing DNS and I want to get k and epsilon after calculation, how should I do?

February 25, 2021, 13:54		#10
ukaya New Member Umut Kaya Join Date: Aug 2020 Posts: 5 Rep Power: 5	It seems like you couldn't copy the function properly to controldict. Also try to run it without all the keywords, but the necessary ones: simpleFoam -postProcess -latestTime