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

Stability problems with kepsilon in external aero

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

Reply
 
LinkBack Thread Tools Display Modes
Old   May 21, 2008, 13:55
Default I'm evaluating OpenFOAM for ex
  #1
New Member
 
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8
edreed is on a distinguished road
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?
edreed is offline   Reply With Quote

Old   May 23, 2008, 04:34
Default Hi Edward, It is a common pra
  #2
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
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
dmoroian is offline   Reply With Quote

Old   May 23, 2008, 11:47
Default I've tried coming up with a ca
  #3
New Member
 
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8
edreed is on a distinguished road
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 is offline   Reply With Quote

Old   May 23, 2008, 14:51
Default With a freestream velocity of
  #4
New Member
 
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8
edreed is on a distinguished road
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?
edreed is offline   Reply With Quote

Old   May 23, 2008, 22:11
Default run checkYplus, what values it
  #5
Senior Member
 
mkraposhin's Avatar
 
Matvej Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 172
Rep Power: 8
mkraposhin is on a distinguished road
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?
mkraposhin is offline   Reply With Quote

Old   May 26, 2008, 01:43
Default Did you used those values (k =
  #6
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
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
dmoroian is offline   Reply With Quote

Old   May 27, 2008, 08:17
Default Average y+ on the surface of t
  #7
New Member
 
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8
edreed is on a distinguished road
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.
edreed is offline   Reply With Quote

Old   May 27, 2008, 09:44
Default Hi Edward, It is surprising t
  #8
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
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
dmoroian is offline   Reply With Quote

Old   May 27, 2008, 10:13
Default Because the definition of char
  #9
New Member
 
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8
edreed is on a distinguished road
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;
}
edreed is offline   Reply With Quote

Old   May 27, 2008, 10:52
Default Ok, I had similar problems wit
  #10
Senior Member
 
dmoroian's Avatar
 
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11
dmoroian is on a distinguished road
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
dmoroian is offline   Reply With Quote

Old   May 27, 2008, 11:00
Default Thanks for the info Dragos, I'
  #11
New Member
 
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8
edreed is on a distinguished road
Thanks for the info Dragos, I'm giving it a shot now.
edreed is offline   Reply With Quote

Old   May 27, 2008, 13:56
Default After further fiddling, I've g
  #12
New Member
 
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8
edreed is on a distinguished road
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.
edreed is offline   Reply With Quote

Old   May 29, 2008, 08:13
Default Hello, Edward. Are you usin
  #13
Member
 
Paulo Alexandre Costa Rocha
Join Date: Mar 2009
Posts: 71
Rep Power: 8
paulo is on a distinguished road
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
paulo is offline   Reply With Quote

Old   May 29, 2008, 08:46
Default Hi Paulo, I am having very
  #14
New Member
 
Ryan Middleton
Join Date: Mar 2009
Posts: 17
Rep Power: 8
ryan_m is on a distinguished road
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
ryan_m is offline   Reply With Quote

Old   May 29, 2008, 08:54
Default I generally do 5 non-orthogona
  #15
New Member
 
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8
edreed is on a distinguished road
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.
edreed is offline   Reply With Quote

Old   May 29, 2008, 08:56
Default Anytime a problem like this ha
  #16
Senior Member
 
Anonymous
Join Date: Mar 2009
Posts: 110
Rep Power: 8
madad2005 is on a distinguished road
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?
madad2005 is offline   Reply With Quote

Old   May 29, 2008, 11:50
Default I'm looking at relatively low
  #17
New Member
 
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8
edreed is on a distinguished road
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.
edreed is offline   Reply With Quote

Old   May 29, 2008, 16:25
Default Hello Ryan and All, Our cyl
  #18
Member
 
Paulo Alexandre Costa Rocha
Join Date: Mar 2009
Posts: 71
Rep Power: 8
paulo is on a distinguished road
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
paulo is offline   Reply With Quote

Old   May 29, 2008, 18:42
Default Hi all, Thanks very much fo
  #19
New Member
 
Ryan Middleton
Join Date: Mar 2009
Posts: 17
Rep Power: 8
ryan_m is on a distinguished road
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
ryan_m is offline   Reply With Quote

Old   May 30, 2008, 10:16
Default Hi Ryan and All, Here we ar
  #20
Member
 
Paulo Alexandre Costa Rocha
Join Date: Mar 2009
Posts: 71
Rep Power: 8
paulo is on a distinguished road
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.
paulo 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
Usefulness of Similarity theory in External Aero Amod Kumar Main CFD Forum 2 January 10, 2007 05:58
boundary conditions for external automotive aero Andrew Berner FLUENT 4 November 2, 2006 12:17
External Aero-boundary condition. Guest FLUENT 0 April 14, 2006 19:31
External Aero Questions Alan FLUENT 4 September 23, 2005 13:58
External aero recommendations leo FLUENT 8 July 1, 2002 02:45


All times are GMT -4. The time now is 05:35.