
[Sponsors] 
September 20, 2007, 03:39 
Hi everyone,
I am strugglin

#1 
Member
Patrick Bourdin
Join Date: Mar 2009
Posts: 40
Rep Power: 8 
Hi everyone,
I am struggling to get simpleFoam working on a purely tetrahedral mesh featuring anisotropic elements in the boundary layer (generated with Gridgen). I tried quite conservative settings with regards to the numerical schemes and the solver parameters (see below), however the computation keeps blowing up (everincreasing continuity error...). I also tried different initializations (internal velocity field = inlet value or converged solution from potentialFoam) but without success. The case at hand is a 3D wing attacked by a flow at an angle of 5 deg. The outer (farfield) boundary has a brick shape: bottom and front faces are inlet (fixedValue velocity and zeroGradient pressure), top and rear faces are pressure outlet (fixedValue pressure and zeroGradient velocity  note that the bc I want there should be velocity outlet, ie zeroGradient for pressure and velocity, but this would make the computation more unstable in the early stages, so defaulted it to pressure outlet). The 2 remaining faces (the side ones) are respectively a symmetry plane (half the wing is modeled) and a slip surface (zero normal component for vectors and zeroGradient for scalars). The turbulence model is SA (inlet value for nuTilda is 5*nu; wall value is zero, Reynolds number based on mean chord is above 300000, y+ of wall adjacent cells should be < 1 based on flatplate correlations). checkMesh was OK, but complained about highly skewed cells though. What else could I tweak in the settings to get the computation converging? Below are the fvSchemes and fvSolutions files I used, as well as a sample of the log file. Cheers, Patrick ddtSchemes { default steadyState; } gradSchemes { default leastSquares; grad(p) leastSquares; grad(U) leastSquares; } divSchemes { default none; div(phi,U) Gauss upwind; div(phi,k) Gauss upwind; 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; } laplacianSchemes { default none; laplacian(nuEff,U) Gauss linear limited 0.333; laplacian((1A(U)),p) Gauss linear limited 0.333; laplacian(DkEff,k) Gauss linear limited 0.333; laplacian(DepsilonEff,epsilon) Gauss linear limited 0.333; laplacian(DREff,R) Gauss linear limited 0.333; laplacian(DnuTildaEff,nuTilda) Gauss linear limited 0.333; } interpolationSchemes { default linear; interpolate(U) linear; } snGradSchemes { default limited 0.333; } fluxRequired { default no; p; } ===================================== solvers { p GAMG { tolerance 1e7; relTol 0.001; smoother GaussSeidel; nPreSweeps 0; nPostSweeps 2; nFinestSweeps 2; cacheAgglomeration true; nCellsInCoarsestLevel 100; agglomerator faceAreaPair; mergeLevels 1; }; U PBiCG { preconditioner DILU; tolerance 1e05; relTol 0.1; }; k PBiCG { preconditioner DILU; tolerance 1e05; relTol 0.1; }; epsilon PBiCG { preconditioner DILU; tolerance 1e05; relTol 0.1; }; R PBiCG { preconditioner DILU; tolerance 1e05; relTol 0.1; }; nuTilda PBiCG { preconditioner DILU; tolerance 1e06; relTol 0.01; }; } SIMPLE { nNonOrthogonalCorrectors 5; pRefCell 0; pRefValue 0; } relaxationFactors { p 0.2; U 0.5; k 0.5; epsilon 0.5; R 0.5; nuTilda 0.5; } ======================================== Time = 27 Lookup gradScheme for grad(U) Lookup divScheme for div((nuEff*dev(grad(U).T()))) Lookup laplacianScheme for laplacian(nuEff,U) Lookup fluxRequired for U Lookup gradScheme for snGradCorr(U) Lookup gradScheme for snGradCorr(U) Lookup gradScheme for snGradCorr(U) Lookup divScheme for div(phi,U) Lookup gradScheme for grad(p) DILUPBiCG: Solving for Ux, Initial residual = 0.0466498, Final residual = 0.00249187, No Iterations 1 DILUPBiCG: Solving for Uy, Initial residual = 0.0354671, Final residual = 0.00175394, No Iterations 1 DILUPBiCG: Solving for Uz, Initial residual = 0.069161, Final residual = 0.0041251, No Iterations 1 Lookup interpolationScheme for interpolate(U) Lookup laplacianScheme for laplacian((1A(U)),p) Lookup fluxRequired for p Lookup gradScheme for snGradCorr(p) GAMG: Solving for p, Initial residual = 3.38756e08, Final residual = 3.38756e08, No Iterations 0 Lookup laplacianScheme for laplacian((1A(U)),p) Lookup fluxRequired for p Lookup gradScheme for snGradCorr(p) GAMG: Solving for p, Initial residual = 3.38756e08, Final residual = 3.38756e08, No Iterations 0 Lookup laplacianScheme for laplacian((1A(U)),p) Lookup fluxRequired for p Lookup gradScheme for snGradCorr(p) GAMG: Solving for p, Initial residual = 3.38756e08, Final residual = 3.38756e08, No Iterations 0 Lookup laplacianScheme for laplacian((1A(U)),p) Lookup fluxRequired for p Lookup gradScheme for snGradCorr(p) GAMG: Solving for p, Initial residual = 3.38756e08, Final residual = 3.38756e08, No Iterations 0 Lookup laplacianScheme for laplacian((1A(U)),p) Lookup fluxRequired for p Lookup gradScheme for snGradCorr(p) GAMG: Solving for p, Initial residual = 3.38756e08, Final residual = 3.38756e08, No Iterations 0 Lookup laplacianScheme for laplacian((1A(U)),p) Lookup fluxRequired for p Lookup gradScheme for snGradCorr(p) GAMG: Solving for p, Initial residual = 3.38756e08, Final residual = 3.38756e08, No Iterations 0 Lookup fluxRequired for p time step continuity errors : sum local = 635.134, global = 0.0084708, cumulative = 0.0724731 Lookup gradScheme for grad(p) Lookup gradScheme for grad(U) Lookup gradScheme for grad(nuTilda) Lookup laplacianScheme for laplacian(DnuTildaEff,nuTilda) Lookup fluxRequired for nuTilda Lookup gradScheme for snGradCorr(nuTilda) Lookup divScheme for div(phi,nuTilda) Lookup ddtScheme for ddt(nuTilda) DILUPBiCG: Solving for nuTilda, Initial residual = 0.14121, Final residual = 0.000208447, No Iterations 3 ExecutionTime = 3756.82 s ClockTime = 3759 s 

September 20, 2007, 15:13 
Does checkMesh report errors o

#2 
Senior Member
BastiL
Join Date: Mar 2009
Posts: 471
Rep Power: 11 
Does checkMesh report errors on your mesh?


September 20, 2007, 15:40 
checkMesh complains about 28 h

#3 
Member
Patrick Bourdin
Join Date: Mar 2009
Posts: 40
Rep Power: 8 
checkMesh complains about 28 highly skewed faces. Is that enough to upset the continuity (it's only 28 in over 7M faces). Perhaps I should give a try to skewLinear wherever linear has been used so far.
Here is the output of checkMesh: /**\  =========    \ / F ield  OpenFOAM: The Open Source CFD Toolbox   \ / O peration  Version: 1.4.1   \ / A nd  Web: http://www.openfoam.org   \/ M anipulation   \**/ Exec : checkMesh . stfwSym_RF00_RA00 Date : Sep 20 2007 Time : 20:17:47 Host : kittyhawk PID : 13253 Root : /users/aexpb/OpenFOAM/aexpb1.4.1/run Case : stfwSym_RF00_RA00 Nprocs : 1 Create time Create polyMesh for time = constant Time = constant Mesh stats points: 661557 edges: 4559608 faces: 7755110 internal faces: 7673122 cells: 3857058 boundary patches: 5 point zones: 0 face zones: 0 cell zones: 0 Number of cells of each type: hexahedra: 0 prisms: 0 wedges: 0 pyramids: 0 tet wedges: 0 tetrahedra: 3857058 polyhedra: 0 Checking topology... Boundary definition OK. Point usage OK. Upper triangular ordering OK. Topological cell zipup check OK. Face vertices OK. Faceface connectivity OK. Number of regions: 1 (OK). Checking patch topology for multiply connected surfaces ... Patch Faces Points Surface wing 56432 28344 ok (not multiply connected) symmPlane 22084 11207 ok (not multiply connected) outflow 1108 602 ok (not multiply connected) inlet 1126 611 ok (not multiply connected) sideWall 1238 658 ok (not multiply connected) Checking geometry... Domain bounding box: (20 0 20) (20 20 20) Boundary openness (3.5074e17 5.11643e16 4.02084e17) OK. Max cell openness = 5.53599e14 OK. Max aspect ratio = 425.218 OK. Minumum face area = 4.03727e11. Maximum face area = 3.65069. Face area magnitudes OK. Min volume = 1.02556e15. Max volume = 2.39606. Total volume = 32000. Cell volumes OK. Mesh nonorthogonality Max: 89.9626 average: 56.9199 *Number of severely nonorthogonal faces: 3053544. Nonorthogonality check OK. <<Writing 3053544 nonorthogonal faces to set nonOrthoFaces Face pyramids OK. ***Max skewness = 10.1112, 28 highly skew faces detected which may impair the quality of the results <<Writing 28 skew faces to set skewFaces *Edges too small, min/max edge length = 3.08586e06 3.42547, number too small: 140120 <<Writing 168479 points on short edges to set shortEdges All angles in faces OK. All face flatness OK. Failed 1 mesh checks. End 

September 20, 2007, 15:51 
Overall, that does not look to

#4 
Senior Member
BastiL
Join Date: Mar 2009
Posts: 471
Rep Power: 11 
Overall, that does not look to bad....
I also had some troubles with simplefoam but finally everything should work better if you try something similar to this: http://www.cfdonline.com/OpenFOAM_D...es/1/4245.html I recommend to try this. This mesh should be possible to run. 

September 20, 2007, 15:53 
In my previous link use the se

#5 
Senior Member
BastiL
Join Date: Mar 2009
Posts: 471
Rep Power: 11 
In my previous link use the settings proposes by Hrvoje Jasak for the car geometry, I missed to mention....


September 21, 2007, 14:12 
Hrv's settings are more a less

#6 
Member
Patrick Bourdin
Join Date: Mar 2009
Posts: 40
Rep Power: 8 
Hrv's settings are more a less the same as mine (see above, first post), mine are expected to be less prone to unstability (stronger limited nonorthogonal correction). Anyways, tried with Hrv's settings, same results: the continuity imbalance wanders in the blue yonder...
Same applies when I do a laminar run (turbulence off in turbulenceProperties) at a much lower Re number. Any other suggestions before I try and do some cosmetic surgery on the mesh (not sure if that will help either)? 

September 21, 2007, 15:10 
Hello Patrick,
I just happe

#7 
Senior Member
Philippose Rajan
Join Date: Mar 2009
Location: Germany
Posts: 530
Rep Power: 16 
Hello Patrick,
I just happened to look through your checkMesh Output, and I am not in the least surprised that your simpleFoam simulation is blowing itself to smithereens! You have 3.85 Million Tetrahedra, and out of these, 3.05 Million cells are severely nonorthogonal (thats almost 80% of your mesh)! You have a maximum nonorthogonality of 89.96 degrees, which is.... lets say.... pretty much as far as it can go :)! And the average nonorthogonality is also quite high at over 56 degrees. I think the nonorthogonal correctors are only meant to improve "mildly" nonorthogonal meshes....not such extreme cases. There also seem to be problems with the range of your edge lengths.... going from xxxe06 to xxxe+0... a huge difference... which could be arising from extreme aspect ratios, and skew. My personal suggestion would be that you schedule a major cosmetic surgery on your mesh as soon as possible :)! Maybe you should first run a repair on your imported geometry before creating the mesh... looks like your geometry has some very small (possibly degenerate) edges. On the other hand.... if you actually do manage to get a working simulation in OpenFOAM with the current mesh that you have, I would be interested to know how you got it working :)! Have a nice weekend! Philippose 

September 21, 2007, 15:13 
Hello Patrick,
Just for en

#8 
Senior Member
Mark Olesen
Join Date: Mar 2009
Location: http://olesenm.github.io/
Posts: 777
Rep Power: 18 
Hello Patrick,
Just for entertainment, have you tried dualizing the mesh and running that? 

September 21, 2007, 15:20 
Have you tried BC's as your s

#9 
New Member
Rajneesh
Join Date: Mar 2009
Posts: 28
Rep Power: 8 
Have you tried BC's as your setup would see in a windtunnel? inlet at the front, outlet at the rear and symmetry on other 4 walls.
Visulaising flowfield just before it blows up may also give you some ideas on where things are going wrong. If you see few elements on your wing surface with high velocity/pressure then problem is likely due to mesh quality. Otherwise changing BC should help. 

September 21, 2007, 15:31 
OK,
The nonorthogonality c

#10 
New Member
Daniel Einstein
Join Date: Mar 2009
Posts: 22
Rep Power: 8 
OK,
The nonorthogonality check says: Mesh nonorthogonality Max: 89.9626 average: 56.9199 *Number of severely nonorthogonal faces: 3053544. Nonorthogonality check OK. Since the check gets a vote of "OK' I am assuming that 90 degreess is maximally orthogonal, in which case, his mesh is pretty good not pretty bad. Am I readingthis backwards, and if so why does the evaluation come out 'OK'? That would be pretty confusing. 

September 23, 2007, 17:42 
Hi guys,
the nonorthogonal

#11 
Member
Patrick Bourdin
Join Date: Mar 2009
Posts: 40
Rep Power: 8 
Hi guys,
the nonorthogonality is so great (almost 90 deg between the cell face normal vector and the displacement vector from owner to neighbour cell centres), because of the meshing strategy: anisotropic tets in the BL and isotropic tets outside (see pics below of the wing section in the symmetry plane). So very close to the wall, where the cell aspect ratio is the greatest, I end up with highly nonorthogonal faces (eg, the hypothenuse shared by 2 triangles in the symmetry plane). However, these socalled anistropic tet elements are supposed to be beneficial in the BL (compared to a squashed isotropic tet), hexas or prisms would be better though. I guess that type of mesh will only work with solvers which obtain the diffusive fluxes at a cell face by interpolating first the stress tensors of the owner and neighbor cells and then carry out the inner product of that interpolated tensor with the face area vector, instead of using a central finite difference (with explicit correction if necessary) to compute directly the normal gradient ("a la" OpenFOAM). Does anyone know how to do that with OpenFOAM (I mean interpolating the entire stress tensor at the cell face instead of using a finite difference scheme to compute the gradient)? 

November 16, 2010, 06:24 

#12 
New Member
Join Date: Jul 2009
Posts: 7
Rep Power: 8 
Hi,
I know this is an old post, but can't find anything more recent covering this topic. I'm having similar issues with Gridgen anisotropic boundary layer meshes diverging. Does anyone know if the OP ever got satisfactory results using anisotets/OpenFOAM? Thanks. 

December 14, 2010, 12:59 

#13 
Senior Member
Travis Carrigan
Join Date: Jul 2010
Location: Arlington, TX
Posts: 122
Rep Power: 6 
Hi,
I managed to put together an aniso case for the NACA 6412 airfoil at 5 degrees angle of attack. I attached airfoilTrex.tar.gz which contains the fvSchemes and fvSolutions dictionaries as well as an image of the residuals demonstrating convergence (probably should be run a little longer). But I thought I would share with you some of the things I did to get this thing running. Mesh details:  30,000 cells  Max aspect ratio=118 in 2D  Initial cell height in BL=0.0001 Simulation details:  2D  NACA 6412  5 degrees AOA  Re=100,000  SpalartAllmaras model with wall functions OpenFOAM details:  version 1.5dev  simpleFOAM Originally I set up the case similar to what was originally posted in this thread. I too had divergence issues. However, I made a few major changes that seemed to make all the difference in the solution and convergence. (take a look at the attached files for more details) Some of the changes: fvSchemes 1. cellLimited leastSquares gradient scheme 2. div(phi,U) divergence scheme set to Gauss linearUpwindV Gauss linear 3. limited 0.5 for both laplacian and snGrad schemes (Note: reconCentral for the divergence and interpolation schemes didn't help at all, it actually caused the solution to diverge) fvSolution 1. I tightened the solver tolerances to 1e8 and relative tolerance to 0 to force the solution to converge to the solver tolerance 2. nNonOrthogonalCorrectors set to 2 Probably the biggest problem I noticed was when the maximum aspect ratio of the aniso cells were well above 100, the solution had trouble converging. It seemed the lower the aspect ratio, hence the higher the overall cell count, the better the convergence. This is pretty significant, but makes some sense. As the aspect ratio of these right angled triangles increase, they become highly nonorthogonal as the angle between the face normal vector and the line connecting adjacent cell centers approaches 90 degrees. However, as the aspect ratio is decreased, this angle becomes smaller and approaches 0 degrees when the cells become isotropic, the ideal case. I found that an aspect ratio of 100 for aniso cells in the boundary layer seems to be a good limit. So in my case my max aspect ratio peaked at about 118. Attachment:  airfoilTrex.tar.gz 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
HOW TO SPECIFY AN ANISOTROPIC SUBSTANCE  VISHNU  CFX  2  December 20, 2007 11:12 
Anisotropic  John  FLUENT  1  June 4, 2007 15:41 
anisotropic thermal conductivity  Lugdi  CDadapco  0  January 15, 2007 09:03 
anisotropic conduction  ramki  Phoenics  2  October 12, 2005 10:36 
Anisotropic Adaptation  Apurva  Main CFD Forum  2  July 10, 2000 17:32 