CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Implement an les turbulence model (http://www.cfd-online.com/Forums/openfoam-programming-development/122711-implement-les-turbulence-model.html)

pante August 26, 2013 09:28

Implement an les turbulence model
 
1 Attachment(s)
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:

Attachment 24851

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::printStack(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::operator/<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 August 26, 2013 09:36

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

chegdan August 26, 2013 09:57

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

Code:

a = b/C
and C could be zero then you can do

Code:

a=b/(C + <some really small number>)
or you can do something like

Code:

a = b/max(SMALL,C)
where SMALL is a very small constant in OpenFOAM. I hope this helps.

pante August 26, 2013 10:36

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 August 27, 2013 13:28

I think you may want to investigate this line further

Code:

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 August 28, 2013 03:15

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 March 11, 2014 06:15

1 Attachment(s)
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 March 17, 2014 12:40

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 August 9, 2014 09:58

Quote:

Originally Posted by pante (Post 479288)
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
Code:

(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 August 9, 2014 23:20

1 Attachment(s)
Hi,pante:
I modified the code to make the expression of Bbeta more easier to read, Attachment 32907
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 August 16, 2014 10:04

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 August 16, 2014 22:42

Quote:

Originally Posted by rajeshkunwar (Post 506248)
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 September 1, 2014 09:38

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 September 1, 2014 12:51

Hi Huang,

this might be helpful:

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

huangxianbei September 1, 2014 20:18

Quote:

Originally Posted by pante (Post 508701)

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 September 1, 2014 22:20

1 Attachment(s)
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
Code:

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
Attachment 33492

huangxianbei September 2, 2014 08:37

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

http://www.cfd-online.com/Forums/mem...ormula-dvm.png

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

Светлана November 24, 2014 00:45

Did you figure it out huangxianbei?

huangxianbei November 24, 2014 02:40

Quote:

Originally Posted by Светлана (Post 520719)
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 December 5, 2014 17:16

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


All times are GMT -4. The time now is 13:04.