[edreed]

I'm evaluating OpenFOAM for external aero applications and have run into a problem while trying to use the k-epsilon turbulence model. I am running the same case in Fluent with k-epsilon no problem, but simpleFoam is immediately unstable.

I am starting from a solution generated by laminar simpleFoam (no turbulence model turned on), so my field velocities are reasonable. I turn on kepsilon and set the boundary conditions on the inlet and the field values to k=e=1 (the same as what I use in Fluent). The first iteration of simpleFoam returns the message:

bounding epsilon, min: -41159.4 max: 1.2112e+06 average: 9541.7

The solution goes wacky from there and simpleFoam crashes within 3 iterations.

Like I said, this exact case works just fine in Fluent using k-epsilon.

Anyone know what I might be doing wrong?

[dmoroian]

Hi Edward,
It is a common practice to set a more appropriate initial value for k and epsilon, even in Fluent. I recommend you to look in the manual for it: k and epsilon initial values. A similar procedure is presented in the Fluent manual (though they suggest a smaller mixing length: l = 0.07*characteristicLength).

I hope this is useful,
Dragos

[edreed]

I've tried coming up with a calculation for initial values of k and epsilon, but it is very unclear to me what should be used for external aero applications. I am not simulating a wind tunnel so I can't just use the inlet turbulence parameters for a specific tunnel.

The main problem I guess is coming up with a characteristic length for an aircraft. I suppose chord of the wing would be appropriate for a fixed wing, but I'm focused on rotocraft fuselages.

[edreed]

With a freestream velocity of 67.909 m/s, assuming 5% of that for U'x=U'y=U'z, I get k = 17.29 m2s-2.

From there assuming 7% of a characteristic length of 1 m, I get epsilon = 158.9 m2s-3.

I apply those as initial values to the field (which already has a reasonable velocity profile from a laminar solution) as well as fixedValue boundary conditions on the inlet.

After 2 iterations, epsilon is bounded and the solution crashes shortly afterward.

Any idea what portion of this I'm screwing up?

[mkraposhin]

run checkYplus, what values it reports?
check boundary file (constant/polyMesh) - what patch type for wall-boundaries? wall?
what initial and boundary conditions are you using?

[dmoroian]

Did you used those values (k = 17.29 m2s-2; epsilon = 158.9 m2s-3) only at the inlet?
You should also initialize the domain with them.

Dragos

[edreed]

Average y+ on the surface of the aircraft is 150ish. I am using wall-boundary patches for the aircraft surface. I have a domain with a flat outlet and a curved inlet (to allow for changing the angle of attack). It looks like a egg with one end cut off flat. The rear patch is set as outlet with p=0 Pa, the curved inlet patch is set as inlet with Ux=67.909 m/s, Uy=Uz=0 m/s. The domain is initialized to Ux=67.909 m/s.

I am setting both k and epsilon at both the inlet patch and the domain.

[dmoroian]

Hi Edward,
It is surprising that you can run the computation without turbulence, and it crashes when you enable it. I would expect an opposite behavior for a turbulent flow (works with turbulence enabled and crashes without).
So my first question is: what Reynolds number do you have, and how do you compute it?
Second, what discretization schemes are you using, especially for the divergence term?

Dragos

[edreed]

Because the definition of characteristic length is somewhat up in the air, I've been using 1 m (the aircraft model is about 3.5 m long). Obviously this is air and I'm just doing standard sea level to start, so nu=1.4607e-5 m2s-1 puts Re=4.65e6.

ddtSchemes
{
default steadyState;
}

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

divSchemes
{
default none;
div(phi,U) Gauss linearUpwind Gauss linear;
div(phi,k) Gauss linearUpwind Gauss linear;
div(phi,epsilon) Gauss linearUpwind Gauss linear;
div(phi,R) Gauss upwind;
div(R) Gauss linear;
div(phi,nuTilda) Gauss upwind;
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;
}

but I've also tried

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

[dmoroian]

Ok, I had similar problems with yours, and the solution was to use

divSchemes
{
<blockquote>default Gauss upwind phi;

}
</blockquote>
and a much better convergence for k and epsilon each iteration:

k PBiCG
{
<blockquote>preconditioner DILU;
tolerance 1e-07;
relTol 0;

};
epsilon PBiCG
{

preconditioner DILU;
tolerance 1e-12;
relTol 0;

};
</blockquote>
My suggestion is to run with Gauss upwind until convergence, and then switch for a higher scheme.

I hope this is helpful,
Dragos

[edreed]

Thanks for the info Dragos, I'm giving it a shot now.

[edreed]

After further fiddling, I've gotten a case to run although epsilon is bounded the entire time.

It wouldn't work with the settings you recommended Dragos, but if I kept the discretization schemes as I had them and upped the initial k and epsilon to 200 and 6000, respectively, I was able to get a solution. Max epsilon in the solution is around 5e7. I got those k and epsilon values from examining a converged Fluent k-e simulation.

I'm starting to think my mesh may be too non-orthogonal to use k-e. Adding additional equations to solve just pushes it over the limit.

Have you had any experience running k-e on polyhedral meshes? I'm doing this case on both a tetrahedral mesh and the polyhedral version of the same mesh, converted using polyDualMesh. The polyhedral version is more unstable and won't run with the settings that allow the tetrahedral case to run.

[paulo]

Hello, Edward.

Are you using non-ortogonal correctors and relaxation?

Here we use them (3 non-orthogonal correctors and a relaxation factor of about 0.3 in the begining of the simulation) in meshes created in Salome (always tetrahedral), and the simulation runs fine.

Good luck!

Paulo A. C. Rocha

[ryan_m]

Hi Paulo,

I am having very similar problems to Edward. I am calculating external aerodynamic flows, and the solution normally runs fine for 100-200 iterations and then epsilon explodes (i turned turbulence off and pressure exploded!) I am using 2 non-orthogonal correctors and under relaxation as 0.7 for all, except pressure with 0.3.

Ive been fiddling around a lot with discretization schemes but not having much luck. Are you able to post your system folder, or give me some hints as to your discretization schemes?

Any help will be much appreciated.

Cheers,
Ryan

[edreed]

I generally do 5 non-orthogonal correctors and I've tried reducing the relaxation on k and epsilon to 0.3 from 0.7 with no luck.

After examining the Fluent solutions and the single OpenFOAM solution I was able to converge using k-epsilon, I think my mesh is just too skewed in places to allow for decent results with a 2-equation model. Maybe I'd have better luck with Spalart-Allmaras, but I think I'll focus on creating better grids as well.

[madad2005]

Anytime a problem like this has occurred for me, mesh quality was normally the culprit. What type of external flow are you aiming for anyway?

[edreed]

I'm looking at relatively low speed (150ish kts) fuselage aerodynamics for rotorcraft. My computational resources limit meshes to generally less than 5-6 million cells and I'm more interested in capturing behavior over a range of pitch/yaw angles (lift curve slope, drag bucket, etc) than I am in answers that absolutely correlate with wind tunnel tests.

[paulo]

Hello Ryan and All,

Our cylinder case is here:

http://www.posmec.ufc.br/~paulo/Open...perc_ke.tar.gz

Please any comments are welcome.

Paulo

[ryan_m]

Hi all,

Thanks very much for all the posts. Paulo, Ill have a good look at your case and see if I can take anything out of it. Looks like Im going to have to add in some more non-orthogonal correctors as a start.

Adriano: you said mesh quality was normally the culprit in your cases. Do you have a recommendation as to the maximum skewness before these problems occur (3D tetrahedral mesh)? I am creating meshes in Gambit and scaling them down in OpenFOAM using transformPoints. My project involves calculating drag coefficients for 3D airfoils and eventually entire aircraft.

Thanks for all the help.

Ryan

[paulo]

Hi Ryan and All,

Here we are also interested in drag for 3D airfoils, and we used the liftdrag utility, as you can see in one of the files in the package.

Unfortunately the results were not good at the start, and I did not have time to investigate yet.

If you have any encouraging results, please let me know.

Sorry for the bad english.

Regards,

Paulo.