[asharma]

Hi Foamers,

I've noticed something unusal in my pisoFoam turbulent simluation, i've been though a couple of other threads regarding relevance of nut and nutilda values, however still not clear in its implementation with the RAS models. I'm currently using the standard k-epsilon model, and the wiki article mentions on how the nut values are modeled. Looking through the code:-

00117     nut_
00118     (
00119         IOobject
00120         (
00121             "nut",
00122             runTime_.timeName(),
00123             mesh_,
00124             IOobject::NO_READ,
00125             IOobject::AUTO_WRITE
00126         ),
00127         autoCreateNut("nut", mesh_)
00128     )
00129 {
00130     nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
00131     nut_.correctBoundaryConditions();
00132 
00133     printCoeffs();
00134 }

It clear how the nut value is being calculated , how can i keep track of how this value is varying in my simulation. My run doesnt seem to be calculating it as opening the nut file at a given time only displays my boundary conditions, i will be requiring nut value to club with my scaler transport model in combination with DT...

thank you in advance,

Regards,
Abhinav

[mturcios777]

I assume that this code is from the turbulence models int he src directory. What you need to do is to add another field in the createFields.H that takes the values of the nut_ field. Creating the field is simple enough:

volScalarField nuTwrite
 (
        IOobject
         (
             "nuTwrite",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        turbulence->nut()
 );

(Assuming you've already created the turbulence model).

You'll probably have to update the field every time step, after the turbulence correction:

turbulence->correct();
nuTwrite = turbulence->nut();

This is all if you want to have nut written out to the time directories of your simulation. If all you need is to access the field, then

turbulence->nut();

should be enough.

[asharma]

Hello Marco,
Thanks for the prompt reply. As suggested i included another field in the createfields.h file, but when i try to compile it, it gives me the following error...

In file included from mypisoFoam.C:46:
createFields.H: In function ‘int main(int, char**)’:
createFields.H:54: error: ‘turbulence’ was not declared in this scope

What am i doing wrong? And for the updation thats to be done in the main solver file right, for example in my case "mypisoFoam.C"...

Regrads,
Abhinav

[boger]

Maybe you should look at something like simpleFoam as an example, which defines "turbulence" as a variable in createFields.H. Note that simpleFoam does not refer to nut explicitly but leaves its management up to the particular turbulence model instead.

[mturcios777]

boger is quite right. Are you basing this custom solver on a previous one, and if so which one. The error you are describing means that the turbulence model has not been created yet, or at all.

alphaEff is a member function of variables of the class turbulenceModel. It may also be that the turbulence model is named differently in your solver. In pisoFoam for example. the declaration is:

autoPtr<incompressible::turbulenceModel> turbulence
(
     incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

,

whereas in the post-processing utility wallHeatFlux it is declared

autoPtr<compressible::RASModel> RASModel
(       
     compressible::RASModel::New
     (   
         rho,
         U,
         phi,
         thermo()
     )
 );

Find out if the turbulence model has been created, and if so what the variable is named (say myTurbulence). Then to get nut, call

myTurbulence->nut()

anytime after that and it should work.

[mturcios777]

And as far as updating the turbulence model, yes that happens in the main solver loop. That way at every time step the entire turbulence model is updated, including nut.

[asharma]

I've based my modified solver on pisofoam, to incorporate scaler transport equation. I had a look through my createfields.H file again and like you said, turbulence is declared as:

autoPtr<incompressible::turbulenceModel> turbulence
    (
        incompressible::turbulenceModel::New(U, phi, laminarTransport)
    );

// part added
volScalarField nuTwrite
    (
        IOobject
         (
             "nuTwrite",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
     turbulence->nut() 
    );

I added the nut term as mentioned earlier. Also in my mail solver 'C' file i update nut values as:-

turbulence->correct();
	nuTwrite = turbulence->nut();

I was able to successfully compile it this time around. I must admit I'm still a little confused on the nuTwrite bit added to the code, as in will my nut values will be stored in nuTwrite folder in my simulation results?

When i tried running my simulation it gave me the following error:

Time = 0.1

Courant Number mean: 0.215709 max: 0.707652
#0 Foam::error:rintStack(Foam::Ostream&) in "/opt/openfoam171/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigSegv::sigSegvHandler(int) in "/opt/openfoam171/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted:
#3 void Foam::add<double, double, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<Foam::typeOfSu m<double, double>::type, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/opt/openfoam171/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#4 Foam::tmp<Foam::GeometricField<Foam::typeOfSum<dou ble, double>::type, Foam::fvPatchField, Foam::volMesh> > Foam:perator+<double, double, Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/opt/openfoam171/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#5 Foam::incompressible::RASModel::nuEff() const in "/opt/openfoam171/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#6 Foam::incompressible::RASModels::kEpsilon::divDevR eff(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&) const in "/opt/openfoam171/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#7
in "/opt/openfoam171/applications/bin/linuxGccDPOpt/mypisoFoam"
#8 __libc_start_main in "/lib/libc.so.6"
#9
in "/opt/openfoam171/applications/bin/linuxGccDPOpt/mypisoFoam"
Segmentation fault

I havn't come across this before, i'm guessing it has something to do with my nut values (or related deceleration) now?

thanks,

Regards,
Abhinav

[asharma]

i re-complied my files but this time i omitted the updation statement to nuTwrite ie:

 turbulence->correct();
	//nuTwrite = turbulence->nut();

my simulation works fine now, but i still don't see how my nut values are being calculated and stored. As Marco mentioned i should be able to get access to the field on adding

myTurbulence->nut()

, what exactly was mean by that? As a test i'm varying my initialized nut values in the 0 directory but it doesn't effect my output velocity and pressure fields?!! In which case can i deduce that its not being used in the simulation run?

Thanks in advance....

Regards,
Abhinav

[mturcios777]

Creating the nutWrite field as an IOobject allows you to have a copy of the fields being written out at every timestep so that you are able to analyze it later in post-processing. In each time directory, you will end up with a nutWrite file that stores all the values of nut for your simulation.

I don't think that the write call is causing the crash, perhaps its symptomatic of it. The stack trace that is printed shows that a segFault occurred when adding who fields (likely nu and nut, since nuEff was called), and following the function calls shows that nuEff was called by divDevRhoeff.

(for more information on how to read traces, check http://openfoamwiki.net/index.php/HowTo_debugging)

On the surface it seems like there is a problem with the turbulence model. kEpsilon based model can have stuff like this happen if the k or epsilon values go wild. Does your solver crash immediately, or does it run for a while before it does?

If it crashed immediately, check the boundary conditions for your turbulence related quantities. If it manages to run for a few timesteps, then reduce your writeInterval (in controlDict) and look at what the turbulence fields (k, epsilon, etc) look like, and if any are going to 0, infinity, or anything else odd (you can do this with paraFoam).

Good luck!

[boger]

Abhinav,

If changes to your nut input file are not effecting your solution, and the nut field is not being written out automatically, then I have to wonder whether the turbulence model is turned on at all. Have you checked the values in the constant/turbulenceProperties and constant/RASproperties files? Does the solver output to the screen/log file show that the turbulence equation(s) are being solved?

[asharma]

nuTwrite = turbulence->nut();

asks the values of nut to be written to the file called nuTwrite right? interestingly the code doesn't crash if i dont have the values written to the file? I've tried this a couple of times but i'm not quite sure whats going on here.
My run crashes and the 2nd time step, so I checked my k and epsilon values, they seemed fine. I did notice that i now have my nuTwrite file writen to my time directory now ( just the one time step) showing me the nut values. I have values shooting up to 16 here!... is this something i should look into?
I also verified the constant/turbulenceProperties

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      turbulenceProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

simulationType  RASModel;

and my constant/RASproperties

RASModel        kEpsilon;

turbulence      on;

printCoeffs     on;

When i start my simulation , the solver displays all the values of the constants

Create time

Create mesh for time = 0

Reading field p

Reading field U

Reading field nut
Reading field T

Reading field source

Reading transportProperties

Reading diffusivity DT

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RASModel
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    sigmaEps        1.3;
}

It does seem to appear that my turbulence is being calculated and the constant values initialized however changing initial nut values doesn't bring about a difference in my velocity pressure output? I'm just worried with my nut values not giving an output, i'm just not confident with my turbulence model.

Thank you Marco and David for your time, your help is much appreciated!

Regards,
Abhinav

[boger]

Frankly, since nut is calculated from k and epsilon, it's always been my feeling that it shouldn't be provided as an input file anyway. Try changing the values of k and epsilon in 0/ -- these should certainly affect your solution.

[mturcios777]

I looked at the tutorial case for pisoFoam, and noticed that it writes nut by default, so there may be something off with your turbulence model that is only being caught when we try copy its nut values to another field.

Besides looking at the initial conditions, there might be something in the way you initialize the turbulence model. Would you be able to post your createFields.H file here for us to look at? I have a hunch about something, but can't be sure until I see.

[asharma]

Here's my createfields.H file on my modified pisofoam model:-

    Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    Info<< "Reading field nut" << endl;
    volScalarField nut
    (
        IOobject
        (
            "nut",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
     );


    Info<< "Reading field T\n" << endl;
    volScalarField T
    (
	IOobject
	(
	    "T",
	    runTime.timeName(),
	    mesh,
	    IOobject::MUST_READ,
	    IOobject::AUTO_WRITE
	    ),
	    mesh
	);
	
     Info<< "Reading field source\n" << endl;
    volScalarField source
    (
	IOobject
	(
	    "source",
	    runTime.timeName(),
	    mesh,
	    IOobject::MUST_READ,
	    IOobject::AUTO_WRITE
	    ),
	    mesh
	);
    
    
    

    Info<< "Reading transportProperties\n" << endl;
    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    );


    Info<< "Reading diffusivity DT\n" << endl;
    dimensionedScalar DT
    (
        transportProperties.lookup("DT")
    );

#   include "createPhi.H"


    label pRefCell = 0;
    scalar pRefValue = 0.0;
    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);


    singlePhaseTransportModel laminarTransport(U, phi);

    autoPtr<incompressible::turbulenceModel> turbulence
    (
        incompressible::turbulenceModel::New(U, phi, laminarTransport)
    );


volScalarField nuTwrite
    (
        IOobject
         (
             "nuTwrite",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
     turbulence->nut() 
    );

Since nut values are being calculated internally through k and epsilon , changing them values should bring about a change in my results. In that case my initial boundary conditions(however differently i initialize them) for nut eventually converge to the same nut values since i've been keeping my k and epsilon values the same. But yes i would like to keep a tab on the values, as i would like to club scaler transport to the pisofoam model by incorporating the nut values in the diffusivity term.

[mturcios777]

Hmmmm, looks pretty standard here. I'm guessing the nut field that is read is similar to the one given in pisoFoam. Are you using the same boundary conditions? Initial values of nut are irrelevant, as nut is calculated from k and epsilon fields. What the 0/nut file does do is specify the boundary conditions. If you look at the cavity tutorial for pisoFoam (which is what I'm guessing where you got your 0/ directory before making modifications) the nut values everywhere are 0, but the boundary conditions show that all the walls are nutWallFunctions.

As one final resort, try changing

nuTWrite = turbulence->nut();

to

turbulence->nut();

If your code still crashes, then there is something wrong with the turbulence model and we need to dig deeper. If your code runs to completion, then there is something wrong with the way nuTWrite is initialized (my fault). I don't know why this is for some fields, but you may need to declare and initialize it with an initial value of 0. You can try change the initialization from

turbulence->nut()

to

mesh,
 dimensionedScalar("nuTWrite",turbulence->nut().dimensions(),0.0)

or something similar.

[asharma]

Yes i've been using the same boundary conditions, except for my inlet which i calculate using groovyBC to account for parabolic profile. Well it worked!! i changed

nuTWrite = turbulence->nut();

to just

turbulence->nut();

Now my run doesn't crash. I do get nuTWrite folder in each time step , but then values are only given for my cell centers at my boundaries, unlike for U or p where i have values for all the cell centers in my domain. However I am satisfied that its calculating nut values which i what i wanted. But shouldn't it display values at all my cell center values?

Thanks a ton Marco and David... really helped me get this up and running!

Regards,
Abhinav

[mturcios777]

No problem. You can comment out the

turbulence->nut()

call now, as it seems that for some reason trying to assign values by the = operator to nuTWrite is what crashes the code. Perhaps it is unnecessary to do this assignment, with it automatically updated by the turbulence model.

Does the nut field that gets written change over time, or does it remain constant? Also, could you post a sample nuTWrite? I'm curious as to the types of fields shown.

[asharma]

Hmm my nut values dont seem to be changing with time.. here's my nuTwrite file 100 sec into my simulation....

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.7.x                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "100";
    object      nuTwrite;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -1 0 0 0 0];

internalField   uniform 16.6734;

boundaryField
{
    canyon
    {
        type            nutWallFunction;
        Cmu             0.09;
        kappa           0.41;
        E               9.8;
        value           nonuniform List<scalar> 
150
(
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.00227466
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.0026869
0.00227466
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
0.00350831
)
;
    }
    inlet
    {
        type            fixedValue;
        value           uniform 0.0001;
    }
    outlet
    {
        type            nutWallFunction;
        Cmu             0.09;
        kappa           0.41;
        E               9.8;
        value           nonuniform List<scalar> 
30
(
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
)
;
    }
    atmosphere
    {
        type            symmetryPlane;
    }
    frontAndback
    {
        type            empty;
    }
}


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

its interesting to note that at higher time steps the last bit of values displayed ie

30
(
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
0.0081488
)
;

isn't shown anymore. Not quite sure whats happening here!.....

Regards,
Abhinav

[mturcios777]

its interesting to note that at higher time steps the last bit of values displayed ie isn't shown anymore. Not quite sure whats happening here!.....

Possibly the values is actually reaching a uniform value to all significant digits. Does the value turn into "uniform", or does it simply stop writing the file?

As far as the constant nut values in time, I'm not sure what to make of it without knowing more about what you are simulating, and even then I may not be able to help. If you have any data to validate against (experimental or other reliable numerical simulation) that would be best. To me it still seems odd that you are specifying nut at the inlet. Is there a physical reason for this?

[asharma]

Yes the value turns into uniform after about 1200 sec into the simulation. I realized there was no point is actually specifying the nut values at the inlet , and as I anyway specify my u,k and epsilon values in terms of a parabolic profile(groovy BC). Nut would then be calculated from these profiles. I will be getting my results validated experimentally, planning on carrying out a monitoring exercise in my city to get data, in the meanwhile I've been looking through a couple of papers to establish my model, including the pollutant dispersion part.