[pbo]

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 (ever-increasing 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 (far-field) 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 flat-plate 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((1|A(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 1e-7;
relTol 0.001;

smoother GaussSeidel;
nPreSweeps 0;
nPostSweeps 2;
nFinestSweeps 2;

cacheAgglomeration true;
nCellsInCoarsestLevel 100;
agglomerator faceAreaPair;
mergeLevels 1;
};
U PBiCG
{
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
};
k PBiCG
{
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
};
epsilon PBiCG
{
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
};
R PBiCG
{
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
};
nuTilda PBiCG
{
preconditioner DILU;
tolerance 1e-06;
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((1|A(U)),p)
Lookup fluxRequired for p
Lookup gradScheme for snGradCorr(p)
GAMG: Solving for p, Initial residual = 3.38756e-08, Final residual = 3.38756e-08, No Iterations 0
Lookup laplacianScheme for laplacian((1|A(U)),p)
Lookup fluxRequired for p
Lookup gradScheme for snGradCorr(p)
GAMG: Solving for p, Initial residual = 3.38756e-08, Final residual = 3.38756e-08, No Iterations 0
Lookup laplacianScheme for laplacian((1|A(U)),p)
Lookup fluxRequired for p
Lookup gradScheme for snGradCorr(p)
GAMG: Solving for p, Initial residual = 3.38756e-08, Final residual = 3.38756e-08, No Iterations 0
Lookup laplacianScheme for laplacian((1|A(U)),p)
Lookup fluxRequired for p
Lookup gradScheme for snGradCorr(p)
GAMG: Solving for p, Initial residual = 3.38756e-08, Final residual = 3.38756e-08, No Iterations 0
Lookup laplacianScheme for laplacian((1|A(U)),p)
Lookup fluxRequired for p
Lookup gradScheme for snGradCorr(p)
GAMG: Solving for p, Initial residual = 3.38756e-08, Final residual = 3.38756e-08, No Iterations 0
Lookup laplacianScheme for laplacian((1|A(U)),p)
Lookup fluxRequired for p
Lookup gradScheme for snGradCorr(p)
GAMG: Solving for p, Initial residual = 3.38756e-08, Final residual = 3.38756e-08, 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

[bastil]

Does checkMesh report errors on your mesh?

[pbo]

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/aexpb-1.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 zip-up check OK.
Face vertices OK.
Face-face 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.5074e-17 5.11643e-16 4.02084e-17) OK.
Max cell openness = 5.53599e-14 OK.
Max aspect ratio = 425.218 OK.
Minumum face area = 4.03727e-11. Maximum face area = 3.65069. Face area magnitudes OK.
Min volume = 1.02556e-15. Max volume = 2.39606. Total volume = 32000. Cell volumes OK.
Mesh non-orthogonality Max: 89.9626 average: 56.9199
*Number of severely non-orthogonal faces: 3053544.
Non-orthogonality check OK.
<<Writing 3053544 non-orthogonal 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.08586e-06 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

[bastil]

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.cfd-online.com/OpenFOAM_D...es/1/4245.html

I recommend to try this. This mesh should be possible to run.

[bastil]

In my previous link use the settings proposes by Hrvoje Jasak for the car geometry, I missed to mention....

[pbo]

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 non-orthogonal 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)?

[philippose]

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 non-orthogonal (thats almost 80% of your mesh)! You have a maximum non-orthogonality of 89.96 degrees, which is.... lets say.... pretty much as far as it can go :-)! And the average non-orthogonality is also quite high at over 56 degrees.

I think the non-orthogonal correctors are only meant to improve "mildly" non-orthogonal meshes....not such extreme cases.

There also seem to be problems with the range of your edge lengths.... going from xxxe-06 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

[olesen]

Hello Patrick,

Just for entertainment, have you tried dualizing the mesh and running that?

[vega]

Have you tried BC's as your setup would see in a wind-tunnel? 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.

[deinstein]

OK,

The non-orthogonality check says:

Mesh non-orthogonality Max: 89.9626 average: 56.9199
*Number of severely non-orthogonal faces: 3053544.
Non-orthogonality 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.

[pbo]

Hi guys,

the non-orthogonality 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 non-orthogonal faces (eg, the hypothenuse shared by 2 triangles in the symmetry plane). However, these so-called 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)?

[Monker1980]

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 aniso-tets/OpenFOAM?

Thanks.

[tcarrigan]

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
- Spalart-Allmaras model with wall functions

OpenFOAM details:
- version 1.5-dev
- 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 1e-8 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 non-orthogonal 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