CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Problems with Turbulence Modeling

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

Reply
 
LinkBack Thread Tools Display Modes
Old   November 26, 2009, 05:48
Default Problems with Turbulence Modeling
  #1
New Member
 
Jignesh
Join Date: Sep 2009
Posts: 6
Rep Power: 8
ezsoal is on a distinguished road
Hi,

I need to simulate flow over a cylinder for the following cases:

1) Inviscid Flow
2) Viscous Subsonic
3) Viscous Subsonic with RAS Turbulence
4) Supersonic

I'm able to simulate the first two (using potentialFoam and icoFoam respectively), but am getting this error while trying to simulate the last two (using pisoFoam and sonicFoam respectiely)

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 RASModel
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
sigmaEps 1.3;
}


Starting time loop

Time = 0.0005

Courant Number mean: 0.000271663 max: 0.00134942
DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 5.73156e-08, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 6.21305e-08, No Iterations 1
DICPCG: Solving for p, Initial residual = 1, Final residual = 0.0913066, No Iterations 31
time step continuity errors : sum local = 8.14153e-05, global = -4.30474e-06, cumulative = -4.30474e-06
DICPCG: Solving for p, Initial residual = 0.0302236, Final residual = 6.3797e-07, No Iterations 49
time step continuity errors : sum local = 7.42487e-08, global = 1.87645e-10, cumulative = -4.30455e-06
#0 Foam::error:rintStack(Foam::Ostream&) in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted:
#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/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 "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/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::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#6 Foam::incompressible::RASModels::kEpsilon::correct () in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#7 main in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/applications/bin/linuxGccDPOpt/pisoFoam"
#8 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#9 _start at /usr/src/packages/BUILD/glibc-2.9/csu/../sysdeps/i386/elf/start.S:122
Floating point exception



Note that the geometry remains the same in all the cases, only the solver and the relevant inputs and boundary types change depending upon the case.

Also, for the supersonic case, the solution blows up (as shown above) after several iteration, whereas for the subsonic turbulent case, it blows up right from the start.

Can anyone please gimme a clue as to what is it that I'm doing wrong?

Regards,
Jignesh
ezsoal is offline   Reply With Quote

Old   November 26, 2009, 08:49
Default
  #2
Member
 
Ulrich Heck
Join Date: Mar 2009
Location: Krefeld, Germany
Posts: 41
Rep Power: 9
ulli is on a distinguished road
Hi Jignesh,

I think there is something wrong with your initial turbulence conditions. This might be caused by your turbulence div-schemes. Waht do you use?

Ulrich
ulli is offline   Reply With Quote

Old   November 26, 2009, 09:44
Default
  #3
New Member
 
Jignesh
Join Date: Sep 2009
Posts: 6
Rep Power: 8
ezsoal is on a distinguished road
Hi Ulrich,

thanks for responding.

I dont have much of a background in CFD. Basically, I've jst copied all the initial data files in the '0' directly (p, U, k, epsilon, nut, nutilda, epsilon) from a similar tutorial and modified them to reflect the patches in my geometry. ( this worked for the previous cases - potentialFoam and icoFoam)

Also, I've simply copied the fvSolution and fvSchemes files from that sample tutorial which solves a flow with RAS Turbulence.

I'm reproducing the contents of these files. I would really appreciate if you could have a look and suggest a way out.

Anticipating response,

Regards,
Jignesh

=====================================
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform 0;

boundaryField
{
down
{
type symmetryPlane;
}

cylinder
{
type symmetryPlane;
}

frontAndBack
{
type empty;
}

fixedWall
{
type zeroGradient;
}

left
{
type zeroGradient;
}

right
{
type fixedValue;
value uniform 0;
}

}

=======================================

FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform (0 0 0);

boundaryField
{
down
{
type symmetryPlane;
}

cylinder
{
type symmetryPlane;
}
frontAndBack

{
type empty;
}

fixedWall
{
type fixedValue;
value uniform (0 0 0);
}

left
{
type fixedValue;
value uniform (5 0 0);
}

right
{
type zeroGradient;
}

}
===============================

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

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

internalField uniform 0.00325;

boundaryField
{
down
{
type symmetryPlane;
}
cylinder
{
type symmetryPlane;
}
frontAndBack
{
type empty;
}
fixedWall
{
type kqRWallFunction;
value uniform 0.00325;
}
left
{
type zeroGradient;
}
right
{
type fixedValue;
value uniform 0;
}
}
==================================

FoamFile
{
version 2.0;
format ascii;
class volSymmTensorField;
object R;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform (0 0 0 0 0 0);

boundaryField
{
down
{
type symmetryPlane;
}

cylinder
{
type symmetryPlane;
}

frontAndBack
{
type empty;
}

fixedWall
{
type kqRWallFunction;
}

left
{
type fixedValue;
value uniform (0 0 0 0 0 0);
}

right
{
type zeroGradient;
}

}
====================================

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

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

internalField uniform 0.000765;

boundaryField
{
down
{
type symmetryPlane;
}

cylinder
{
type symmetryPlane;
}

frontAndBack
{
type empty;
}

fixedWall
{
type epsilonWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 0.000765;
}

right
{
type zeroGradient;
}

left
{
type fixedValue;
value uniform 14.855;
}
}

=====================================

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

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

internalField uniform 0;

boundaryField
{
down
{
type symmetryPlane;
}

cylinder
{
type symmetryPlane;
}

frontAndBack
{
type empty;
}

fixedWall
{
type nutWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 0;
}

left
{
type calculated;
value uniform 0;
}

right
{
type calculated;
value uniform 0;
}
}
=========================================

FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object nuTilda;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform 0;

boundaryField
{
down
{
type symmetryPlane;
}

cylinder
{
type symmetryPlane;
}

frontAndBack
{
type empty;
}

fixedWall
{
type kqRWallFunction;
}

left
{
type fixedValue;
value uniform 0;
}

right
{
type zeroGradient;
}

}
===============================================

FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default Euler;
}

gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss limitedLinearV 1;
div(phi,k) Gauss limitedLinear 1;
div(phi,epsilon) Gauss limitedLinear 1;
div(phi,R) Gauss limitedLinear 1;
div(R) Gauss linear;
div(phi,nuTilda) Gauss limitedLinear 1;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
default none;
laplacian(nuEff,U) Gauss linear corrected;
laplacian((1|A(U)),p) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian(DREff,R) Gauss linear corrected;
laplacian(DnuTildaEff,nuTilda) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
interpolate(U) linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p ;
}
==========================
ezsoal is offline   Reply With Quote

Old   November 26, 2009, 10:46
Default
  #4
Member
 
Ulrich Heck
Join Date: Mar 2009
Location: Krefeld, Germany
Posts: 41
Rep Power: 9
ulli is on a distinguished road
Hi Jignesh

I think you should not force k to zero at your patch "right". Maybe you should try zero gradient instead.

If this doesn't help. Consider more stable schemes for turbulence, e.g.
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;

and switch to more accurate schemes later.
But in most cases, when the solution crashes so soon, there is something wrong with the initial conditions or bcs. You can also set turbulence Intensities and mixingLenght for inlets
e.g. k:

type turbulentIntensityKineticEnergyInlet;
intensity 0.01;
value $internalField;

e.g. omega:
type turbulentMixingLengthFrequencyInlet;
mixingLength 0.04;
k k;
value $internalField;

For epsilon is turbulentMixingLengthDissipationInlet I think.

Hope this helps

Ulli
ulli is offline   Reply With Quote

Old   November 26, 2009, 16:12
Default
  #5
New Member
 
Jignesh
Join Date: Sep 2009
Posts: 6
Rep Power: 8
ezsoal is on a distinguished road
Hey Ulrich,

Thanks for the suggestion; it worked! i got the simulation running flawlessly by NOT forcing k to be equal to zero! now I jst need to ensure that the results that I have are indeed correct...

however, the problem with the supersonic case seems to be something different, since it is blowing up suddenly after 0.0027secs:

Time = 0.00273952

Courant Number mean: 0.00110348 max: 0.0129647
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Ux, Initial residual = 0.970821, Final residual = 4.52189e-06, No Iterations 27
DILUPBiCG: Solving for Uy, Initial residual = 0.950914, Final residual = 8.50896e-06, No Iterations 26
DILUPBiCG: Solving for e, Initial residual = 0.956971, Final residual = 5.481e-06, No Iterations 27
DILUPBiCG: Solving for p, Initial residual = 0.100512, Final residual = 7.64695e-09, No Iterations 13
DILUPBiCG: Solving for p, Initial residual = 0.0364326, Final residual = 9.77965e-09, No Iterations 7
DILUPBiCG: Solving for p, Initial residual = 0.00127526, Final residual = 2.52387e-09, No Iterations 9
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 2.60677e-09, global = 1.43791e-09, cumulative = 2.09749e-08
DILUPBiCG: Solving for p, Initial residual = 0.0276749, Final residual = 2.05278e-10, No Iterations 15
DILUPBiCG: Solving for p, Initial residual = 0.019402, Final residual = 1.06515e-09, No Iterations 7
DILUPBiCG: Solving for p, Initial residual = 0.00143388, Final residual = 7.04595e-09, No Iterations 5
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 8.05887e-08, global = -4.09812e-08, cumulative = -2.00063e-08
DILUPBiCG: Solving for epsilon, Initial residual = 0.415028, Final residual = 6.03136e-09, No Iterations 35
bounding epsilon, min: -4.23037e+26 max: 6.89055e+26 average: 2.43139e+23
DILUPBiCG: Solving for k, Initial residual = 0.980292, Final residual = 2.48219e-09, No Iterations 34
bounding k, min: -3.11137e+16 max: 2.40184e+17 average: 7.80979e+15
ExecutionTime = 6421.04 s ClockTime = 6739 s

Time = 0.00273956

Courant Number mean: 0.108861 max: 259.099
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG: Solving for Ux, Initial residual = 0.949675, Final residual = 8.48593e-06, No Iterations 26
DILUPBiCG: Solving for Uy, Initial residual = 0.902717, Final residual = 8.30647e-06, No Iterations 26
DILUPBiCG: Solving for e, Initial residual = 0.954222, Final residual = 6.66186e-06, No Iterations 26


Maximum number of iterations exceeded#0 Foam::error:rintStack(Foam::Ostream&) in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::error::abort() in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Foam::ePsiThermo<Foam:ureMixture<Foam::constTran sport<Foam::specieThermo<Foam::eConstThermo<Foam:: perfectGas> > > > >::calculate() in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libbasicThermophysicalModels.so"
#3 Foam::ePsiThermo<Foam:ureMixture<Foam::constTran sport<Foam::specieThermo<Foam::eConstThermo<Foam:: perfectGas> > > > >::correct() in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/lib/linuxGccDPOpt/libbasicThermophysicalModels.so"
#4 main in "/home/jignesh/OpenFOAM/OpenFOAM-1.6/applications/bin/linuxGccDPOpt/sonicFoam"
#5 __libc_start_main in "/lib/tls/i686/cmov/libc.so.6"
#6 _start at /usr/src/packages/BUILD/glibc-2.9/csu/../sysdeps/i386/elf/start.S:122


From function specieThermo<thermo>::T(scalar f, scalar T0, scalar (specieThermo<thermo>::*F)(const scalar) const, scalar (specieThermo<thermo>::*dFdT)(const scalar) const) const
in file /home/dm2/henry/OpenFOAM/OpenFOAM-1.6/src/thermophysicalModels/specie/lnInclude/specieThermoI.H at line 68.

FOAM aborting

Aborted


The Courant No is SUDDENLY shooting up by a factor of 20,000 in jst one iteration. Infact the results till t=0.0027 seem to be fairly decent when viewed in paraFoam. but I want to run the simulation till 0.01 to ensure it has reached steady-state; and I have no clue as to WHAT is causing this blowup. Any ideas about trouble shooting? Has it got anything to do with the thermophysical model?

i'm reproducing the files for the supersonic case to make things more clear:

Again, many thanks for the help!

Regards,
Jignesh


========================
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform 100000;

boundaryField
{
inlet
{
type fixedValue;
value uniform 100000;
}

outlet
{
type waveTransmissive;
field p;
phi phi;
rho rho;
psi psi;
gamma 1.3;
fieldInf 100000;
lInf 1;
value uniform 100000;
}

topwall
{
type zeroGradient;
}

symplane
{
type symmetryPlane;
}

cylinder
{
type zeroGradient;
}

frontAndBack
{
type empty;
}
}
================================

FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField uniform (650 0 0);

boundaryField

{
topwall
{
type supersonicFreestream;
pInf 100000;
TInf 300;
UInf (650 0 0);
gamma 1.28418;
value uniform (650 0 0);
}

symplane
{
type symmetryPlane;
}

frontAndBack
{
type empty;
}

cylinder
{
type fixedValue;
value uniform (0 0 0);
}

inlet
{
type fixedValue;
value uniform (650 0 0);
}

outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (0 0 0);
}

}

================================

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

dimensions [0 0 0 1 0 0 0];

internalField uniform 300;

boundaryField

{
inlet
{
type fixedValue;
value uniform 300;
}

outlet
{
type inletOutlet;
inletValue uniform 300;
value uniform 300;
}

symplane
{
type symmetryPlane;
}

topwall
{
type inletOutlet;
inletValue uniform 300;
value uniform 300;
}

cylinder
{
type zeroGradient;
}

frontAndBack
{
type empty;
}
}
===============================

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

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

internalField uniform 1000;

boundaryField
{
inlet
{
type fixedValue;
value uniform 1000;
}
outlet
{
type inletOutlet;
inletValue uniform 1000;
value uniform 1000;
}
topwall
{
type inletOutlet;
inletValue uniform 1000;
value uniform 1000;
}
symplane
{
type symmetryPlane;
}
cylinder
{
type compressible::kqRWallFunction;
value uniform 1000;
}
frontAndBack
{
type empty;
}
}
===================================

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

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

internalField uniform 0;

boundaryField
{
inlet
{
type calculated;
value uniform 0;
}
outlet
{
type calculated;
value uniform 0;
}
topwall
{
type calculated;
value uniform 0;
}
symplane
{
type symmetryPlane;
}
cylinder
{
type mutWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 0;
}
frontAndBack
{
type empty;
}
}

========================================

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

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

internalField uniform 25000;

boundaryField
{
inlet
{
type fixedValue;
value uniform 25000;
}
outlet
{
type inletOutlet;
inletValue uniform 25000;
value uniform 25000;
}
topwall
{
type inletOutlet;
inletValue uniform 25000;
value uniform 25000;
}
symplane
{
type symmetryPlane;
}
cylinder
{
type compressible::epsilonWallFunction;
Cmu 0.09;
kappa 0.41;
E 9.8;
value uniform 25000;
}
frontAndBack
{
type empty;
}
}

=================================================

FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
default Euler;
}

gradSchemes
{
default Gauss linear;
}

divSchemes
{
default none;
div(phi,U) Gauss limitedLinearV 1;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(phi,R) Gauss upwind;
div(R) Gauss linear;
div(phid,p) Gauss limitedLinear 1;
div(phiU,p) Gauss limitedLinear 1;
div(phi,e) Gauss limitedLinear 1;
div((muEff*dev2(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
default none;
laplacian(muEff,U) Gauss linear corrected;
laplacian(DkEff,k) Gauss linear corrected;
laplacian(DREff,R) Gauss linear corrected;
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
laplacian((rho*(1|A(U))),p) Gauss linear corrected;
laplacian(alphaEff,e) Gauss linear corrected;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default corrected;
}

fluxRequired
{
default no;
p;
}
ezsoal 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
Discussion: Reason of Turbulence!! Wen Long Main CFD Forum 3 May 15, 2009 09:52
Y plus problems with 3D modeling Míša Milakovová FLUENT 0 April 1, 2008 12:09
turbulence modeling error at a stagnation point erdem NUMECA 2 August 15, 2006 15:40
steady turbulence modeling kerem Main CFD Forum 0 April 24, 2006 17:04
CFD - Trends and Perspectives Jonas Larsson Main CFD Forum 16 August 7, 1998 16:27


All times are GMT -4. The time now is 16:02.