CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Implement an les turbulence model

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 1 Post By chegdan
  • 1 Post By cfdonline2mohsen

Reply
 
LinkBack Thread Tools Display Modes
Old   August 26, 2013, 09:28
Default Implement an les turbulence model
  #1
New Member
 
Join Date: Feb 2013
Posts: 8
Rep Power: 4
pante is on a distinguished road
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!
pante is offline   Reply With Quote

Old   August 26, 2013, 09:36
Default
  #2
Senior Member
 
Bernhard
Join Date: Sep 2009
Location: Delft
Posts: 790
Rep Power: 12
Bernhard is on a distinguished road
Check if you are not dividing by zero somewhere, maybe due to initialization.
Bernhard is offline   Reply With Quote

Old   August 26, 2013, 09:57
Default
  #3
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 546
Rep Power: 18
chegdan will become famous soon enough
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.
songwukong likes this.
__________________
Dan

Find me on twitter @dancombest and LinkedIn
chegdan is offline   Reply With Quote

Old   August 26, 2013, 10:36
Default
  #4
New Member
 
Join Date: Feb 2013
Posts: 8
Rep Power: 4
pante is on a distinguished road
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..
pante is offline   Reply With Quote

Old   August 27, 2013, 13:28
Default
  #5
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 546
Rep Power: 18
chegdan will become famous soon enough
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.
__________________
Dan

Find me on twitter @dancombest and LinkedIn
chegdan is offline   Reply With Quote

Old   August 28, 2013, 03:15
Default
  #6
New Member
 
Join Date: Feb 2013
Posts: 8
Rep Power: 4
pante is on a distinguished road
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 is offline   Reply With Quote

Old   March 11, 2014, 06:15
Default
  #7
New Member
 
Join Date: Feb 2013
Posts: 8
Rep Power: 4
pante is on a distinguished road
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.
Attached Files
File Type: zip Vreman.zip (6.9 KB, 42 views)
pante is offline   Reply With Quote

Old   March 17, 2014, 12:40
Default
  #8
Member
 
Florian Ries
Join Date: Feb 2014
Location: Darmstadt, Germany
Posts: 88
Rep Power: 3
itchy is on a distinguished road
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
itchy is offline   Reply With Quote

Old   August 9, 2014, 09:58
Default
  #9
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 266
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by pante View Post
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

Last edited by huangxianbei; August 9, 2014 at 23:23.
huangxianbei is offline   Reply With Quote

Old   August 9, 2014, 23:20
Default
  #10
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 266
Rep Power: 4
huangxianbei is on a distinguished road
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
huangxianbei is offline   Reply With Quote

Old   August 16, 2014, 10:04
Default
  #11
New Member
 
Rajesh Kumar
Join Date: Apr 2009
Posts: 23
Rep Power: 8
rajeshkunwar is on a distinguished road
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.
rajeshkunwar is offline   Reply With Quote

Old   August 16, 2014, 22:42
Default
  #12
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 266
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by rajeshkunwar View Post
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

Last edited by huangxianbei; August 19, 2014 at 21:42.
huangxianbei is offline   Reply With Quote

Old   September 1, 2014, 09:38
Default
  #13
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 266
Rep Power: 4
huangxianbei is on a distinguished road
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
Is is possible to use OFstream in a turbulence model?
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
huangxianbei is offline   Reply With Quote

Old   September 1, 2014, 12:51
Default
  #14
New Member
 
Join Date: Feb 2013
Posts: 8
Rep Power: 4
pante is on a distinguished road
Hi Huang,

this might be helpful:

Improved implementation of dynamic Smagorinsky
pante is offline   Reply With Quote

Old   September 1, 2014, 20:18
Default
  #15
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 266
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by pante View Post
Hi Huang,

this might be helpful:

Improved implementation of dynamic Smagorinsky
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 is offline   Reply With Quote

Old   September 1, 2014, 22:20
Default
  #16
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 266
Rep Power: 4
huangxianbei is on a distinguished road
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
DVM-new.tar.gz

Last edited by huangxianbei; September 2, 2014 at 08:30.
huangxianbei is offline   Reply With Quote

Old   September 2, 2014, 08:37
Default
  #17
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 266
Rep Power: 4
huangxianbei is on a distinguished road
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
huangxianbei is offline   Reply With Quote

Old   November 24, 2014, 00:45
Default
  #18
Member
 
Join Date: Oct 2013
Posts: 30
Rep Power: 3
Светлана is on a distinguished road
Did you figure it out huangxianbei?
Светлана is offline   Reply With Quote

Old   November 24, 2014, 02:40
Default
  #19
Senior Member
 
Huang Xianbei
Join Date: Sep 2013
Location: CAU,China
Posts: 266
Rep Power: 4
huangxianbei is on a distinguished road
Quote:
Originally Posted by Светлана View Post
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
huangxianbei is offline   Reply With Quote

Old   December 5, 2014, 17:16
Smile
  #20
Senior Member
 
cfdonline2mohsen's Avatar
 
Mohsen KiaMansouri
Join Date: Jan 2010
Location: CFD Lab
Posts: 118
Rep Power: 7
cfdonline2mohsen is on a distinguished road
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
PonchO likes this.
__________________
“If you have an apple and I have an apple and we exchange these apples then you and I will still each have one apple. But if you have an idea and I have an idea and we exchange these ideas, then each of us will have two ideas.”
cfdonline2mohsen is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
gamma-ReTheta turbulence model for predicting transitional flows FelixL OpenFOAM Programming & Development 77 June 20, 2015 10:20
How to decide to Turbulence model shipman OpenFOAM 2 August 18, 2013 03:00
Wrong calculation of nut in the kOmegaSST turbulence model FelixL OpenFOAM Bugs 27 March 27, 2012 09:02
Turbulence Model and limitation to Reynolds number qascapri FLUENT 0 January 24, 2011 11:48
y+ for LES turbulence model Pawan FLUENT 0 December 7, 2007 01:40


All times are GMT -4. The time now is 07:51.