[shogan50]

I'm playing with a wedge geometry which is meshed with tets. (netgen 1D..3D) in Salome. At the sharp corner of the wedge at the outlet, potential foam generates very large magnitude velocities. I'm beginning to suspect it is due to the tet mesh. I've played with various BC's. The phenomenon seems to occur at both inlet and outlet if the velocity isn't completely constrained.
Suggestions?
My BC's:

inlet:
p:zeroGradient
U0 0 -100);
p: 1e5
U:
outlet
{
type inletOutlet;
value uniform (0 0 0);
inletValue uniform ( 0 0 0);

}

[shogan50]

I modeled this as a half pipe with symmetryPlane. The outlet has the same phenomenon - U values on the order of 70000 in all directions. This partially rules out my previous theory - that it has to do with the point in the wedge. Suggestions?

[sega]

I made bad experience with tet-meshes, but they are not uncommon.

If nothing is wrong with your boundary condition itself you could try to increase the number of non-orthogonal-correctors.
This is the entry in fvSolution called nNonOrthogonalCorrectors.

Increasing this entry to much may result in a severe increase of computational time!
But with the right value you will be able to overcome the problems concerning correct flux calculations when dealing with tet-elements.

Unfortunately there is no guidline how to choose the right value. I suppose you will have to try ...

[shogan50]

Thanks very much for the help. Apparently this falls under the SIMPLE heading in the fvSolution file. It was set to zero and increasing it decreases the phenomenon. At 3 it starts to look good. 10 seems to approach an error asymptote. Below is an output snippet.

I was also confused by functionality in paraview. I had tried this previously. Paraview wasn't automatically changing the scaling, so I didn't notice the change.




Calculating potential flow
DICPCG: Solving for p, Initial residual = 1, Final residual = 0.000951584, No Iterations 61
DICPCG: Solving for p, Initial residual = 0.909937, Final residual = 0.000886357, No Iterations 43
DICPCG: Solving for p, Initial residual = 0.775395, Final residual = 0.000722846, No Iterations 54
DICPCG: Solving for p, Initial residual = 0.847058, Final residual = 0.000837037, No Iterations 39
DICPCG: Solving for p, Initial residual = 0.758512, Final residual = 0.000746753, No Iterations 53
DICPCG: Solving for p, Initial residual = 0.897108, Final residual = 0.000881784, No Iterations 34
DICPCG: Solving for p, Initial residual = 0.89887, Final residual = 0.000770348, No Iterations 55
DICPCG: Solving for p, Initial residual = 0.870719, Final residual = 0.000765984, No Iterations 32
DICPCG: Solving for p, Initial residual = 0.638767, Final residual = 0.000635883, No Iterations 53
DICPCG: Solving for p, Initial residual = 0.26727, Final residual = 0.000256407, No Iterations 32
DICPCG: Solving for p, Initial residual = 0.080843, Final residual = 7.82863e-05, No Iterations 51
continuity error = 0.0420761
Interpolated U error = 1.70791e-07
ExecutionTime = 1.62 s ClockTime = 1 s

[alberto]

Hello,

did you check the skewness of the mesh?

Your pressure equation has a hard time to converge, and after 11 correctors the initial residual is still around 0.1, which means "not converged" ;-)

Best,

[shogan50]

I'm not sure if this is the exact mesh (may have made some changes), but potential foam shows similar results. Here is the checkMesh output. Is there a good description of skewness somewhere? I'm using netgen 1d-2d-3d in Salome to mesh a step output solidworks model.

scotth@scotth-ubuntu:~/OpenFOAM/scotth-1.6/run/trumpetLES$ checkMesh
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.6                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 1.6-53b7f692aa41
Exec   : checkMesh
Date   : Apr 03 2010
Time   : 09:34:13
Host   : scotth-ubuntu
PID    : 10279
Case   : /home/scotth/OpenFOAM/scotth-1.6/run/trumpetLES
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           6907
    faces:            65092
    internal faces:   58628
    cells:            30930
    boundary patches: 5
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     0
    prisms:        0
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    30930
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch               Faces    Points   Surface topology                  
    wall                3415     1804     ok (non-closed singly connected)  
    inlet               624      350      ok (non-closed singly connected)  
    outlet              365      210      ok (non-closed singly connected)  
    exh                 53       38       ok (non-closed singly connected)  
    symmetry            2007     1072     ok (non-closed singly connected)  

Checking geometry...
    Overall domain bounding box (-0.109217 0 0) (0.109217 0.109151 0.3048)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (5.4753e-19 -5.58808e-20 -4.65674e-19) OK.
    Max cell openness = 1.3833e-16 OK.
    Max aspect ratio = 4.69247 OK.
    Minumum face area = 1.40048e-05. Maximum face area = 0.000151843.  Face area magnitudes OK.
    Min volume = 2.60542e-08. Max volume = 5.24061e-07.  Total volume = 0.00433954.  Cell volumes OK.
    Mesh non-orthogonality Max: 51.3832 average: 16.7662
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.678317 OK.

Mesh OK.

[alberto]

FYI, similar problems on meshes with tets were described here

http://www.cfd-online.com/Forums/ope...tml#post252778

Best,

[makaveli_lcf]

Regarding my ivestigations at http://www.cfd-online.com/Forums/openfoam-solving/73744-overestimated-temperature-values.html:

I made calculations with foamCalc:

1. foamCalc div phi maximal error of order 10^-4
2. foamCalc div U maximal errors of order 10^2!!!

For U I set "div(U) Gauss linear;" in fvSchemes
Does that a reason for early described errors?

[alberto]

Hello,

you should look at div(phi) since phi is the conservative flux, and yes, that error is telling you the solution is not converging well.

Try reducing the tolerances on the pressure linear solver to 10^-15 or 10^-20 (it will slow down the solution a lot, usually with these strict requirements conjugate gradients methods are faster than GAMG in my experience). In addition, consider more non-orthogonal correctors (the number depends on how your specific convergence behaviour).

Best,
Alberto

[makaveli_lcf]

Thanks, I'm trying your guides right now...
Should div(U) also go to zero or is that something apart from flux conservation?

[alberto]

It should be a lot smaller than what you have. The difference between div(U) and div(phi) is due to the interpolation. The flux phi is the variable responsible of ensuring the conservation, while you can look at U as a secondary variable.

Best,

[makaveli_lcf]

Ok! Trying to treat it more, because I found that most errors for div(U) are at the same cells, which appear as an numerical heat sources in my case...

[alberto]

Lovely skewed cells... :-D

[makaveli_lcf]

Some one could please explain me following result?
I got this picture in icoFoam cavity case, tolerance for U and p is 1e-30, steady state achieved.

divU.jpg

Picture represents logarithmic plot of the interpolated velocities divergence.
Utility, which I used to achieve such result is also attached.

divCheck.tar.gz

You can compile utility with wmake command and use then:

divCheck -time YourTime


Waiting for your comments!

[alberto]

If you compute fvc::div(phi), which is what you have to check, since phi is the conservative flux, you obtain (cavity tutorial, time = 0.5, probably not at perfect steady state):

Divergens absolute max/min/avg values: 0.0004 0 4.729e-05

Best,

[makaveli_lcf]

Alberto,

thank you for your comment! I also have checked it before with "foamCalc div phi" and calculated the flow for 5 seconds, sorry I didn't mention that.
But my question is - why flux, interpolated from U cell values, does not correspond to continuity equation

div U = 0?

Is this a PISO / SIMPLE feature? Or perhaps some other scheme apart from "linear" should be used for interpolate(U)? And do you mean, that if "div phi = 0" I should not worry about mass conservation?

[alberto]

Hi,

in OpenFOAM you enforce the mass conservation by imposing

div(phi) = 0

where phi is the mass or volumetric flux (meaning the velocity interpolated on faces, dotted with the surface + terms due to co-located grids). This is what you actually do also on codes based on the staggered grid arrangement, but there you do not risk any confusion, because the velocity does not need to be interpolated, since you have it already on the face.

From the previous equation, you compute a pressure that ensures the the flux phi is conservative after the correction is applied, and then, you use phi in the linearization of the momentum equation and in the evaluation of the other convective terms (for example in scalar transport equations).

So, strictly speaking, when you check if the mass conservation is enforced correctly, you should check div(phi). See for example continuityErrs.H

I would say that, if you want to calculate the difference between div(phi), with phi coming from the pEqn, and the interpolated U, you should compute

div(fvc::interpolate(U) & mesh.Sf())

after the velocity has been corrected with the updated pressure field. At this point you can calculate the difference of the correct flux and the interpolated one.

Best,

[makaveli_lcf]

Thank you, Alberto! Now it is more clear for me.