CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   Stability problems with kepsilon in external aero (

edreed May 21, 2008 14:55

I'm evaluating OpenFOAM for ex
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 May 23, 2008 05:34

Hi Edward, It is a common pra
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,

edreed May 23, 2008 12:47

I've tried coming up with a ca
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 May 23, 2008 15:51

With a freestream velocity of
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 May 23, 2008 23:11

run checkYplus, what values it
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 May 26, 2008 02:43

Did you used those values (k =
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.


edreed May 27, 2008 09:17

Average y+ on the surface of t
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 May 27, 2008 10:44

Hi Edward, It is surprising t
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?


edreed May 27, 2008 11:13

Because the definition of char
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.

default steadyState;

default Gauss linear;
grad(p) Gauss linear;
grad(U) Gauss linear;

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;

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;

default linear;
interpolate(U) linear;

default corrected;

default no;

but I've also tried

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 May 27, 2008 11:52

Ok, I had similar problems wit
Ok, I had similar problems with yours, and the solution was to use
<blockquote>default Gauss upwind phi;

and a much better convergence for k and epsilon each iteration:
<blockquote>preconditioner DILU;
tolerance 1e-07;
relTol 0;
epsilon PBiCG
preconditioner DILU;
tolerance 1e-12;
relTol 0;

My suggestion is to run with Gauss upwind until convergence, and then switch for a higher scheme.

I hope this is helpful,

edreed May 27, 2008 12:00

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

edreed May 27, 2008 14:56

After further fiddling, I've g
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 May 29, 2008 09:13

Hello, Edward. Are you usin
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 May 29, 2008 09:46

Hi Paulo, I am having very
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.


edreed May 29, 2008 09:54

I generally do 5 non-orthogona
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 May 29, 2008 09:56

Anytime a problem like this ha
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 May 29, 2008 12:50

I'm looking at relatively low
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 May 29, 2008 17:25

Hello Ryan and All, Our cyl
Hello Ryan and All,

Our cylinder case is here:

Please any comments are welcome.


ryan_m May 29, 2008 19:42

Hi all, Thanks very much fo
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.


paulo May 30, 2008 11:16

Hi Ryan and All, Here we ar
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.



All times are GMT -4. The time now is 06:48.