
[Sponsors] 
May 21, 2008, 13:55 
I'm evaluating OpenFOAM for ex

#1 
New Member
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8 
I'm evaluating OpenFOAM for external aero applications and have run into a problem while trying to use the kepsilon turbulence model. I am running the same case in Fluent with kepsilon 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 kepsilon. Anyone know what I might be doing wrong? 

May 23, 2008, 04:34 
Hi Edward,
It is a common pra

#2 
Senior Member
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11 
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 

May 23, 2008, 11:47 
I've tried coming up with a ca

#3 
New Member
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8 
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. 

May 23, 2008, 14:51 
With a freestream velocity of

#4 
New Member
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8 
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 m2s2.
From there assuming 7% of a characteristic length of 1 m, I get epsilon = 158.9 m2s3. 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? 

May 23, 2008, 22:11 
run checkYplus, what values it

#5 
Senior Member
Matvey Kraposhin
Join Date: Mar 2009
Location: Moscow, Russian Federation
Posts: 225
Rep Power: 10 
run checkYplus, what values it reports?
check boundary file (constant/polyMesh)  what patch type for wallboundaries? wall? what initial and boundary conditions are you using?
__________________
Winter OpenFOAM days in Moscow  http://opencloudconf.ru/en 

May 26, 2008, 01:43 
Did you used those values (k =

#6 
Senior Member
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11 
Did you used those values (k = 17.29 m2s2; epsilon = 158.9 m2s3) only at the inlet?
You should also initialize the domain with them. Dragos 

May 27, 2008, 08:17 
Average y+ on the surface of t

#7 
New Member
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8 
Average y+ on the surface of the aircraft is 150ish. I am using wallboundary 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. 

May 27, 2008, 09:44 
Hi Edward,
It is surprising t

#8 
Senior Member
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11 
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 

May 27, 2008, 10:13 
Because the definition of char

#9 
New Member
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8 
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.4607e5 m2s1 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((1A(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; } 

May 27, 2008, 10:52 
Ok, I had similar problems wit

#10 
Senior Member
Dragos
Join Date: Mar 2009
Posts: 647
Rep Power: 11 
Ok, I had similar problems with yours, and the solution was to use
divSchemes} </blockquote> and a much better convergence for k and epsilon each iteration: k PBiCG}; epsilon PBiCG { preconditioner DILU;}; </blockquote> My suggestion is to run with Gauss upwind until convergence, and then switch for a higher scheme. I hope this is helpful, Dragos 

May 27, 2008, 11:00 
Thanks for the info Dragos, I'

#11 
New Member
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8 
Thanks for the info Dragos, I'm giving it a shot now.


May 27, 2008, 13:56 
After further fiddling, I've g

#12 
New Member
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8 
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 ke simulation. I'm starting to think my mesh may be too nonorthogonal to use ke. Adding additional equations to solve just pushes it over the limit. Have you had any experience running ke 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. 

May 29, 2008, 08:13 
Hello, Edward.
Are you usin

#13 
Member
Paulo Alexandre Costa Rocha
Join Date: Mar 2009
Posts: 71
Rep Power: 8 
Hello, Edward.
Are you using nonortogonal correctors and relaxation? Here we use them (3 nonorthogonal 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 

May 29, 2008, 08:46 
Hi Paulo,
I am having very

#14 
New Member
Ryan Middleton
Join Date: Mar 2009
Posts: 17
Rep Power: 8 
Hi Paulo,
I am having very similar problems to Edward. I am calculating external aerodynamic flows, and the solution normally runs fine for 100200 iterations and then epsilon explodes (i turned turbulence off and pressure exploded!) I am using 2 nonorthogonal 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 

May 29, 2008, 08:54 
I generally do 5 nonorthogona

#15 
New Member
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8 
I generally do 5 nonorthogonal 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 kepsilon, I think my mesh is just too skewed in places to allow for decent results with a 2equation model. Maybe I'd have better luck with SpalartAllmaras, but I think I'll focus on creating better grids as well. 

May 29, 2008, 08:56 
Anytime a problem like this ha

#16 
Senior Member
Anonymous
Join Date: Mar 2009
Posts: 110
Rep Power: 8 
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?


May 29, 2008, 11:50 
I'm looking at relatively low

#17 
New Member
Edward Reed
Join Date: Mar 2009
Posts: 12
Rep Power: 8 
I'm looking at relatively low speed (150ish kts) fuselage aerodynamics for rotorcraft. My computational resources limit meshes to generally less than 56 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.


May 29, 2008, 16:25 
Hello Ryan and All,
Our cyl

#18 
Member
Paulo Alexandre Costa Rocha
Join Date: Mar 2009
Posts: 71
Rep Power: 8 
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 

May 29, 2008, 18:42 
Hi all,
Thanks very much fo

#19 
New Member
Ryan Middleton
Join Date: Mar 2009
Posts: 17
Rep Power: 8 
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 nonorthogonal 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 

May 30, 2008, 10:16 
Hi Ryan and All,
Here we ar

#20 
Member
Paulo Alexandre Costa Rocha
Join Date: Mar 2009
Posts: 71
Rep Power: 8 
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. 

Thread Tools  
Display Modes  


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 Aeroboundary 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 