[hylleman]

Hi FOAM:ers,

I'm trying to get nice lift/drag prediction on 2D airfoils at low AoA:s (0-5 degrees) and a Re of 10^6, using the simpleFoam solver. I've used some low-Re turbulence models, especially SST k-omega, and some different meshes, but am consequently getting to high drag, an overprediction of 100-150%. From the threads here on this subject I understand that this is normal for fully turbulent modells, and that a transitional model may be needed.

Thus I have now turned to the kklomega model, which I believe is the only transitional model implemented in OpenFOAM. The problem is that I can't even get this model to run. In the first iteration, the solver ends with a floating point exception. I would very much appriciate help on making this case run.

These are my BC:s,

Kt:
inlet fixedValue 3.4e-4
outlet inletOutlet inletValue 3.4e-4
airfoil fixedValue 1e-11

Kl:
inlet fixedValue 1e-11
outlet inletOutlet inletValue 1e-11
airfoil fixedValue 1e-11

omega:
inlet fixedValue 0.25
outlet inletOutlet 0.25
airfoil omegaWallFunction

nut:
inlet/outlet calculated
airfoil nutLowReWallFunction

p:
inlet zeroGradient
outlet fixedValue 0
airfoil zeroGradient

U:
inlet fixedValue (15.11 0 0)
outlet inletOutlet inletValue (15.11 0 0)
airfoil fixedValue (0 0 0)

These are my numerical schemes,

divSchemes
{
default none;
div(phi,U) bounded Gauss limitedLinearV 1;
div(phi,kt) bounded Gauss limitedLinear 1;
div(phi,kl) bounded Gauss limitedLinear 1;
div(phi,omega) bounded Gauss limitedLinear 1;
div((nuEff*dev(T(grad(U))))) Gauss linear;
}


I've been using this wonderful mesher by Mads here on the forum,
http://hvirvel.dk/airfoilmesher/
Should probably also mention that I'm using OpenFOAM 2.3.0. I would very much appriciate help on making this case run.

Regards,
Karl

EDIT: After reading the article a bit closer, I changed omega to zeroGradient at the airfoil. Not the problem though.

EDIT2: The error message might be helpful. It suggests a division by zero in the kklomega::correct function, if I understand it correctly. I cannot see what might be zero at the start though...

smoothSolver:  Solving for Ux, Initial residual = 0.020463, Final residual = 0.00138025, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.0726345, Final residual = 0.00486526, No Iterations 1
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.00940194, No Iterations 2
time step continuity errors : sum local = 0.00295125, global = -1.54808e-06, cumulative = -1.54808e-06
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/usr/lib/libc.so.6"
#3  Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#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&) at ??:?
#5  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#6  Foam::incompressible::RASModels::kkLOmega::correct() at ??:?
#7  
 in "/opt/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/bin/simpleFoam"
#8  __libc_start_main in "/usr/lib/libc.so.6"
#9  
 in "/opt/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/bin/simpleFoam"

[hylleman]

I managed to get this working, after a lot of fiddeling. Unfortunatly, I can't say exactly what the problem was since I changed a lot (schemes, bcs, relaxation factors...). Now I get really nice force prediction on the airfoil, so the fiddeling was worth it

[MichaelD]

Hi Karl,

I'm also interested in trying the kklomega model, and I've run into the same problem. Could you post the setup you used to get it working?

Thanks in advance!

[tareqkh]

Originally Posted by hylleman

I managed to get this working, after a lot of fiddeling. Unfortunatly, I can't say exactly what the problem was since I changed a lot (schemes, bcs, relaxation factors...). Now I get really nice force prediction on the airfoil, so the fiddeling was worth it

Hey,

Could you please provide us more explications since I am working on the same problem?

[papercut]

I get the same instant error

smoothSolver:  Solving for Ux, Initial residual = 0.020463, Final residual = 0.00138025, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 0.0726345, Final residual = 0.00486526, No Iterations 1
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.00940194, No Iterations 2
time step continuity errors : sum local = 0.00295125, global = -1.54808e-06, cumulative = -1.54808e-06
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/usr/lib/libc.so.6"
#3  Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
#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&) at ??:?
#5  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&, Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#6  Foam::incompressible::RASModels::kkLOmega::correct() at ??:?
#7  
 in "/opt/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/bin/simpleFoam"
#8  __libc_start_main in "/usr/lib/libc.so.6"
#9  
 in "/opt/OpenFOAM/OpenFOAM-2.3.0/platforms/linux64GccDPOpt/bin/simpleFoam"

Anyone knows how to deal with this?

[alexeym]

Hi,

If you take a look at correct method of kkLOmega, programming there always is rather defensive, like this:

const volScalarField lambdaT(sqrt(kT)/(omega_ + omegaMin_));

to avoid division be zero. Except several places, here's one of them:

    const volScalarField nutl
    (
        min
        (
            C11_*fTaul(lambdaEff, ktL)*omega*sqr(lambdaEff)
          * sqrt(ktL)*lambdaEff/nu()
          + C12_*BetaTS(ReOmega)*ReOmega*sqr(y_)*omega
        ,
            0.5*(kl_ + ktL)/sqrt(S2)
        )
    );

where S2 is

const volScalarField S2(2.0*magSqr(symm(gradU)));

Why your S2 goes to zero somewhere? I don't know, can you show your initial and boundary conditions?

Update. Though if you have got diverging solution and something that should be positive has become negative, then all these omega_ + omegaMin_ can become zero.

[malv83]

After 8 years, there is a new version (or new model) of the k-kl-omega model.

There are a few problems with the k-kl-omega model in the farfield. One of them is the growth of Laminar Kinetic energy when separation occurs. Lopez and Walters have a paper (have not been published yet) correcting this issue:

Maurin Lopez. D. K. Walters. “A recommended correction to the k-kl-omega transition sensitive eddy-viscosity model”. Journal of Fluid Engineering.

This correction has to be made to the 2008 k-kl-omega model from now on.

Now, Lopez and Walters also developed a new transitional model (k-omega-v2) as an alternative to the k-kl-omega one. This new model has more capabilities (it is more reliable) than the k-kl-omega model, especially in the farfield computations. Fortunately the paper for this new model is already publish.

Maurin Lopez. D. K. Walters. “Prediction of transitional and fully turbulent free shear flows using an alternative to the laminar kinetic energy approach”. Journal of Turbulence, Vol 17, Iss. 3, 2016.

If you see the papers, you will immediately see how the k-kl-omega model is not good for free shear flows, and how the new model corrects all those issues. From now on, k-kl-omega users have to start using the new k-omega-v2 model.