[pante]

Dear Foamers,

I' m trying to implement an eddy-viscosity subgrid-scale model for turbulent shear flow, proposed by Vreman (http://www.vremanresearch.nl/Vreman-...bgridmodel.pdf).

I started from the Smagorinsky model and changed it as is shown in the files attached below:

Vreman.zip

I try to run pitzDaily tutorial with pisoFoam after compiling the turbulence model with wmake libso (no errors in compilation) with vreman turbulence model (I added libs ("libVreman.so"); in system/controlDict file) and i get the following:

pisoFoam
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.2.x |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 2.2.x-1ef4d132e582
Exec : pisoFoam
Date : Aug 26 2013
Time : 16:24:35
Host : "dromis.polytech.ucy.ac.cy"
PID : 28772
Case : /home/pantelis/vreman
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

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

Create mesh for time = 0

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type LESModel
Selecting LES turbulence model Vreman
Selecting LES delta type cubeRootVol
--> FOAM Warning :
From function cubeRootVolDelta::calcDelta()
in file cubeRootVolDelta/cubeRootVolDelta.C at line 52
Case is 2D, LES is not strictly applicable

#0 Foam::error:rintStack(Foam::Ostream&) in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2
at sigaction.c:0
#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#4 void Foam::divide<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libincompressibleRASModels.so"
#5 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam:perator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libincompressibleRASModels.so"
#6 Foam::incompressible::LESModels::Vreman::updateSub GridScaleFields(Foam::GeometricField<Foam::Tensor< double>, Foam::fvPatchField, Foam::volMesh> const&) in "/home/pantelis/OpenFOAM/pantelis-2.2.x/platforms/linux64GccDPOpt/lib/libVreman.so"
#7 Foam::incompressible::LESModels::Vreman::Vreman(Fo am::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::transportModel&, Foam::word const&, Foam::word const&) in "/home/pantelis/OpenFOAM/pantelis-2.2.x/platforms/linux64GccDPOpt/lib/libVreman.so"
#8 Foam::incompressible::LESModel::adddictionaryConst ructorToTable<Foam::incompressible::LESModels::Vre man>::New(Foam::GeometricField<Foam::Vector<double >, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::transportModel&, Foam::word const&) in "/home/pantelis/OpenFOAM/pantelis-2.2.x/platforms/linux64GccDPOpt/lib/libVreman.so"
#9 Foam::incompressible::LESModel::New(Foam::Geometri cField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::transportModel&, Foam::word const&) in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libincompressibleLESModels.so"
#10 Foam::incompressible::turbulenceModel::addturbulen ceModelConstructorToTable<Foam::incompressible::LE SModel>::NewturbulenceModel(Foam::GeometricField<F oam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::transportModel&, Foam::word const&) in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libincompressibleLESModels.so"
#11 Foam::incompressible::turbulenceModel::New(Foam::G eometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::transportModel&, Foam::word const&) in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/lib/libincompressibleTurbulenceModel.so"
#12
in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/pisoFoam"
#13 __libc_start_main in "/lib64/libc.so.6"
#14
in "/share/apps/centFOAM/OpenFOAM/OpenFOAM-2.2.x/platforms/linux64GccDPOpt/bin/pisoFoam"
Floating point exception


Can anybody explain me why i get this error?

thanks!

[Bernhard]

Check if you are not dividing by zero somewhere, maybe due to initialization.

[chegdan]

Floating Point Exceptions are almost always associated with a divide by zero. Looking at your implementation, if you have a gradU of zero then you will run in to some issues. I'm away from my Linux machine to correct your code but if you have a relation like

a = b/C

and C could be zero then you can do

a=b/(C + <some really small number>)

or you can do something like

a = b/max(SMALL,C)

where SMALL is a very small constant in OpenFOAM. I hope this helps.

[pante]

Thanks Bernhard and Daniel for your replies..
I initialize the fields with a converged solution (using smagorinsky turbulence model) and then I compute mag(Grad(U)) (the denominator in the equation for nuSgs) in the whole domain. I found the min value of this 6.74. So the denominator is not zero..

[chegdan]

I think you may want to investigate this line further

 nuSgs_ = 2.5*sqrt(2*ck_/ce_)*ck_*sqr(delta())*sqrt(Bbeta(gradU)/a(gradU));

As, if one comments out this line the error will disappear, albeit give the wrong answer of course. Looking further, ce_ is defined in another header as a non-zero constant and if you print it out to the screen it will be a positive non-zero value (so this isn't it). Bbeta and a may not be initializes properly as Bernhard suggested and since Foam::divide shows up in you error message....you may want to look further. good luck.

[pante]

ok guys, i get it work!

I correct:

nuSgs_ = 2.5*sqrt(2*ck_/ce_)*ck_*sqrt(mag(Bbeta(gradU)))/max(a(gradU),small1_);

where

small1_=1e-10

Thanks again!

[pante]

Hi again guys..
Although the results i get when i use vreman model for a pipe flow are quite accurate, the speed of the computations is very slow: ~ 32 s per time step for a mesh of 700k cells on 8 processors (the respective time when i use dynamic smagorinsky model is ~ 8s). I suppose my code need improvements and optimization. Can anybody check my code and suggest improvements?
I attach the files below.

[itchy]

Hi Pante,

which dynamic smagorinsky model do you use?? homogeneousDynSmagorinsky or one with local averaging??

I also want to implement the model. At the moment I check the paper you posted and after that I will check your code.

A first thing:
if you use the max-function, pherhaps this makes your code slow (Its like a clipping). If you want to have it faster try to add small1_ in the counter and dump the max(...,small1_).

kind regards
Florian

[huangxianbei]

Originally Posted by pante

Hi again guys..
Although the results i get when i use vreman model for a pipe flow are quite accurate, the speed of the computations is very slow: ~ 32 s per time step for a mesh of 700k cells on 8 processors (the respective time when i use dynamic smagorinsky model is ~ 8s). I suppose my code need improvements and optimization. Can anybody check my code and suggest improvements?
I attach the files below.

Hi,pante:
I looked into your code and found that the calculation of Bbeta is complex

(pow(delta(),4)*0.5*(sqr(tr(T(gradU)&gradU))-tr((T(gradU)&gradU)&(T(gradU)&gradU))));

I think it should be simplified, however, I haven't got further idea . By the way, is it possible to use tensor components in the code ? say beta.component and so on?

Xianbei

[huangxianbei]

Hi,pante:
I modified the code to make the expression of Bbeta more easier to read, forum.tar.gz
however, the speed doesn't change at all!
here is a comparison of time require in one time steep loop(centrifugal pump)
Vreman-my 324.91s converged in 19 steps
Vreman-Pante 324.43s converged in 19 steps
Smagorinsky 280.88s converged in 20 steps

Xianbei

[rajeshkunwar]

Hello FOAMERS,

As per the formulation of the Vreman (2004), if mag(gradU) is zero, Sgs Viscosity is consistently zero. Therefore, by introducing small1_ in the code, you are increasing the the SGS viscosity.

[huangxianbei]

Originally Posted by rajeshkunwar

Hello FOAMERS,

As per the formulation of the Vreman (2004), if mag(gradU) is zero, Sgs Viscosity is consistently zero. Therefore, by introducing small1_ in the code, you are increasing the the SGS viscosity.

Hi,Rajesh:
The small is used to prevent 0 in the Denominator， if not, the calculation may be terminated by this kind of error.

PS:Now I understand your formula of Bbeta and find it's the best way to implement this model,especially when dynamic model is considered

Xianbei

[huangxianbei]

Hi,Pante:
I get the idea from your Vreman model's code and want to implement the dynamic Vreman model proposed by You(2006). However, the output nuSgs file only contain the values on the boundaries. May I draw your attention to here
http://www.cfd-online.com/Forums/ope...nce-model.html
The dynamic model I implemented is in the 3th thread. Could you please take a look at and give me some advice?

Thank you very much.

Xianbei

[pante]

Hi Huang,

this might be helpful:

http://www.cfd-online.com/Forums/ope...agorinsky.html

[huangxianbei]

Originally Posted by pante

Hi Huang,

this might be helpful:

http://www.cfd-online.com/Forums/ope...agorinsky.html

Hi,Pante:
Thank you. In fact, I'm using that version as in the link. I also look into that dynamic Smagorinsky code, however, I still can't figure out what's wrong in my code. The main problem is that nuSgs can't output as a field. However, the calculation seems not be affected, it can get similar result as the dynamic Smagorinsky model.This makes the problem more harder to fix Any comments will be appreciated.

Xianbei

[huangxianbei]

Hi,Pante:
I rewrite the code based on dynamic Smagorinsky and the Vreman model you implemented. However, things seems not change at all. Here is the output nuSgs in a processor

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

internalField   uniform -0;

boundaryField
{
    in
    {
        type            calculated;
        value           nonuniform 0();
    }
    out
    {
        type            calculated;
        value           uniform -0;
    }
    shroud
    {
        type            nuSgsUSpaldingWallFunction;
        kappa           0.41;
        E               9.8;
        value           uniform 0;
    }
    hub
    {
        type            nuSgsUSpaldingWallFunction;
        kappa           0.41;
        E               9.8;
        value           uniform 0;
    }
    blade
    {
        type            zeroGradient;
    }
    peri1
    {
        type            cyclic;
    }
    peri2
    {
        type            cyclic;
    }
    procBoundary1to0
    {
        type            processor;
        value           uniform -0;
    }
    procBoundary1to0throughperi1
    {
        type            processorCyclic;
        value           uniform -0;
    }
}

And here is the new version of the mdoel
DVM-new.tar.gz

[huangxianbei]

As in the new version, I put all the formulas to the .C file to make it more easy to read and use without go back to .H file.

Here is the formulas used in dynamic Vreman model



As in the code, I use cD to represent Cv*II
Thus, cD is a volScalarField. aa is magSqr(alpha)
and 'fil' corresponding to the filtered quantities.
I adopt your formula about Bbeta and adapt it to the filtered Bbeta.

I think the formulas in the code is quite right

while,still and still

I can't output the right nuSgs field

Xianbei

[Svetlana]

Did you figure it out huangxianbei?

[huangxianbei]

Originally Posted by Светлана

Did you figure it out huangxianbei?

Not clear yet, but you can output full nuSgs by bounding the Cv>0 instead of (nuSgs+nu)>0

[cfdonline2mohsen]

Dear Huang Xianbei

Please take a look at the following links at the Rostock university which have implemented the following SGS models:

• dynamic mixed SGS model (Zhang, Vreman) with extension to scalar transport

• dynamic mixed version of Meneveau's lagrangian turbulence model with extension to scalar transport

They might be helpful.


http://www.lemos.uni-rostock.de/en/d...nfoam-content/


https://github.com/LEMOS-Rostock/LEM...mpressible/LES