[ezsoal]

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

[ulli]

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

[ezsoal]

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 ;
}
==========================

[ulli]

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

[ezsoal]

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;
}