[MartinBlx]

Hi everybody,

I'm facing a weird problem and I haven't managed to fix it on my own so I'm asking your help.

My geometry is basically two chambers connected by a helical channel, the ratio between the two cross section is roughly 100.
My inlet is located at the bottom of the upper chamber, which means the velocity of my fluid (water) is multiplied by 100 in the channel.

I have carried simpleFoam (and pisoFoam) simulations on different meshes.

simpleFoam converged. However, when analyzing the flow at the outlet, it is obvious that the solution is not good:

https://drive.google.com/file/d/0B_T...ew?usp=sharing

https://drive.google.com/file/d/0B_T...ew?usp=sharing

My mesh:

https://drive.google.com/file/d/0B_T...ew?usp=sharing

the checkMesh:

Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
    points:           57486
    faces:            155293
    internal faces:   137764
    cells:            49544
    faces per cell:   5.9150856
    boundary patches: 11
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     46348
    prisms:        1761
    wedges:        0
    pyramids:      0
    tet wedges:    223
    tetrahedra:    0
    polyhedra:     1212
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            4   788
            5   424

Checking topology...
    Boundary definition OK.
    Cell to face addressing 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
                    maxY        0        0                        ok (empty)
                    minX        0        0                        ok (empty)
                    maxX        0        0                        ok (empty)
                    minY        0        0                        ok (empty)
                    minZ        0        0                        ok (empty)
                    maxZ        0        0                        ok (empty)
                   Inlet     2784     2840  ok (non-closed singly connected)
               ChamberIn     4243     4401  ok (non-closed singly connected)
                 Channel     3413     4080  ok (non-closed singly connected)
              ChamberOut     4321     4427  ok (non-closed singly connected)
                  Outlet     2768     2840  ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (-0.047994251 -0.047981246 -4.2513276e-08) (0.047955685 0.04799059 0.024506538)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (2.1704086e-17 -5.0874776e-17 8.0606393e-16) OK.
    Max cell openness = 3.3265621e-16 OK.
    Max aspect ratio = 5.0262167 OK.
    Minimum face area = 1.217387e-07. Maximum face area = 5.1785612e-06.  Face area magnitudes OK.
    Min volume = 2.6254044e-10. Max volume = 5.7125854e-09.  Total volume = 0.00011854518.  Cell volumes OK.
    Mesh non-orthogonality Max: 62.319649 average: 5.4554404
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 2.232373 OK.
    Coupled point location match (average 0) OK.

Mesh OK.

End

The fvSchemes:

ddtSchemes
{
    default             steadyState;
}

gradSchemes
{
    default             Gauss linear;
}

divSchemes
{
    default             none;
    div(phi,U)          bounded Gauss upwind;
    div((nuEff*dev(T(grad(U)))))    Gauss linear;
}

laplacianSchemes
{
    default             Gauss linear limited corrected 0.333;
}

interpolationSchemes
{
    default             linear;
}

snGradSchemes
{
    default             limited corrected 0.333;
}

fluxRequired
{
    default             no;
    p;
}

The fvSolution:

solvers
{
    p
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.1;
        smoother         GaussSeidel;
        nPreSweeps       0;
        nPostSweeps      2;
        cacheAgglomeration on;
        agglomerator     faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels      1;
    }

    U
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-8;
        relTol           0.1;
        nSweeps          1;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;

    residualControl
    {
        p               1e-3;
        U               1e-4;
    }
}

relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
    }
}

cache
{
    grad(U);
}

The log file for the last iterations (it converged):

Time = 985

smoothSolver:  Solving for Ux, Initial residual = 6.38073e-05, Final residual = 5.95454e-06, No Iterations 5
smoothSolver:  Solving for Uy, Initial residual = 4.56813e-05, Final residual = 4.02243e-06, No Iterations 5
smoothSolver:  Solving for Uz, Initial residual = 0.000115959, Final residual = 9.89525e-06, No Iterations 5
GAMG:  Solving for p, Initial residual = 5.72742e-06, Final residual = 1.15765e-07, No Iterations 2
time step continuity errors : sum local = 0.00562212, global = 0.000201768, cumulative = -1.96636
ExecutionTime = 70.03 s  ClockTime = 72 s

Time = 986

smoothSolver:  Solving for Ux, Initial residual = 8.78342e-05, Final residual = 7.97117e-06, No Iterations 5
smoothSolver:  Solving for Uy, Initial residual = 5.06628e-05, Final residual = 4.44554e-06, No Iterations 5
smoothSolver:  Solving for Uz, Initial residual = 0.000128446, Final residual = 1.02145e-05, No Iterations 5
GAMG:  Solving for p, Initial residual = 5.30479e-06, Final residual = 1.93503e-07, No Iterations 2
time step continuity errors : sum local = 0.00939721, global = 0.000386174, cumulative = -1.96597
ExecutionTime = 70.08 s  ClockTime = 73 s

Time = 987

smoothSolver:  Solving for Ux, Initial residual = 6.53228e-05, Final residual = 6.50554e-06, No Iterations 5
smoothSolver:  Solving for Uy, Initial residual = 4.70409e-05, Final residual = 3.94105e-06, No Iterations 5
smoothSolver:  Solving for Uz, Initial residual = 9.76812e-05, Final residual = 7.99976e-06, No Iterations 5
GAMG:  Solving for p, Initial residual = 4.73726e-06, Final residual = 1.7753e-07, No Iterations 2
time step continuity errors : sum local = 0.00861289, global = 0.000316105, cumulative = -1.96566
ExecutionTime = 70.14 s  ClockTime = 73 s


SIMPLE solution converged in 987 iterations

End

Finalising parallel run

I don't know what is the origin of the problem because I successfully carried similar simulations. Is that a meshing problem? or a problem with the schemes???

I'm a bit confused.

[Gerrit]

Hi,

have you tried to use a finer mesh? As you have propably high Re numbers, you should think about the right turbulence model for your case.
Based on that it might make sense to mesh finer at the walls and in front of the outlet.

Hope that helps

[ssss]

Increase the nNonOrthogonalCorrectors and reduce the residual control because your values are too high. Moreover are you sure that your flow is not turbulent?

[MartinBlx]

Here a new mesh I generated (roughly 6 000 000 cells) where I refined in the previous divergence area:

https://drive.google.com/file/d/0B_T...ew?usp=sharing

but the checkMesh gives some bad results (5 failed checks):

Mesh stats
    points:           6171314
    faces:            17650093
    internal faces:   17175464
    cells:            5753314
    faces per cell:   6.05313
    boundary patches: 5
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     5418331
    prisms:        52953
    wedges:        0
    pyramids:      0
    tet wedges:    364
    tetrahedra:    0
    polyhedra:     281666
    Breakdown of polyhedra by number of faces:
        faces   number of cells
            4   27801
            5   13491
            6   7905
            7   124544
            8   45159
            9   55767
           10   19
           11   4
           12   5678
           15   1152
           18   146

Checking topology...
    Boundary definition OK.
    Cell to face addressing 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
                   Inlet    34141    35141  ok (non-closed singly connected)
               ChamberIn   122844   123527  ok (non-closed singly connected)
                 Channel   112119   136825  ok (non-closed singly connected)
              ChamberOut   122845   123515  ok (non-closed singly connected)
                  Outlet    82680    83122  ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (-0.0489987 -0.049 -1.8272e-07) (0.049 0.049 0.0245002)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-2.00751e-16 1.09013e-16 4.71144e-14) OK.
 ***High aspect ratio cells found, Max aspect ratio: 6.63345e+192, number of cells 12632
  <<Writing 12632 cells with high aspect ratio to set highAspectRatioCells
    Minimum face area = 2.60377e-11. Maximum face area = 4.33051e-07.  Face area magnitudes OK.
 ***Zero or negative cell volume detected.  Minimum negative volume: -1.2117e-11, Number of negative volume cells: 12632
  <<Writing 12632 zero volume cells to set zeroVolumeCells
    Mesh non-orthogonality Max: 180 average: 10.0422
   *Number of severely non-orthogonal (> 70 degrees) faces: 16484.
 ***Number of non-orthogonality errors: 40913.
  <<Writing 57397 non-orthogonal faces to set nonOrthoFaces
 ***Error in face pyramids: 95678 faces are incorrectly oriented.
  <<Writing 65404 faces with incorrect orientation to set wrongOrientedFaces
 ***Max skewness = 102532, 5329 highly skew faces detected which may impair the quality of the results
  <<Writing 5329 skew faces to set skewFaces
    Coupled point location match (average 0) OK.

Failed 5 mesh checks.

End

I still tried to run simpleFoam and it crashed before the end of the first iteration with the following error message:

[0] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[0] #1  Foam::sigFpe::sigHandler(int) at ??:?
[0] #2   in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #3  Foam::GAMGSolver::scale(Foam::Field<double>&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const at ??:?
[0] #4  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
[0] #5  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[0] #6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
[0] #7  Foam::fvMatrix<double>::solve(Foam::dictionary const&) at ??:?
[0] #8  
[0]  at ??:?
[0] #9  
[0]  at ??:?
[0] #10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[0] #11  
[0]  at ??:?
[mp13305-All-Series:10557] *** Process received signal ***
[mp13305-All-Series:10557] Signal: Floating point exception (8)
[mp13305-All-Series:10557] Signal code:  (-6)
[mp13305-All-Series:10557] Failing at address: 0x3e90000293d
[mp13305-All-Series:10557] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x36d40) [0x7f8f2110bd40]
[mp13305-All-Series:10557] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x39) [0x7f8f2110bcc9]
[mp13305-All-Series:10557] [ 2] /lib/x86_64-linux-gnu/libc.so.6(+0x36d40) [0x7f8f2110bd40]
[mp13305-All-Series:10557] [ 3] /opt/openfoam231/platforms/linux64GccDPOpt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver5scaleERNS_5FieldIdEES3_RKNS_9lduMatrixERKNS_10FieldFieldIS1_dEERKNS_8UPtrListIKNS_17lduInterfaceFieldEEERKS2_h+0xbe) [0x7f8f22256e9e]
[mp13305-All-Series:10557] [ 4] /opt/openfoam231/platforms/linux64GccDPOpt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver6VcycleERKNS_7PtrListINS_9lduMatrix8smootherEEERNS_5FieldIdEERKS8_S9_S9_S9_S9_S9_RNS1_IS8_EESD_h+0x222d) [0x7f8f2225b74d]
[mp13305-All-Series:10557] [ 5] /opt/openfoam231/platforms/linux64GccDPOpt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver5solveERNS_5FieldIdEERKS2_h+0x4ae) [0x7f8f2225c32e]
[mp13305-All-Series:10557] [ 6] /opt/openfoam231/platforms/linux64GccDPOpt/lib/libfiniteVolume.so(_ZN4Foam8fvMatrixIdE15solveSegregatedERKNS_10dictionaryE+0x132) [0x7f8f23fa32c2]
[mp13305-All-Series:10557] [ 7] simpleFoam(_ZN4Foam8fvMatrixIdE5solveERKNS_10dictionaryE+0x14f) [0x44678f]
[mp13305-All-Series:10557] [ 8] simpleFoam() [0x446aba]
[mp13305-All-Series:10557] [ 9] simpleFoam() [0x41b9a1]
[mp13305-All-Series:10557] [10] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5) [0x7f8f210f6ec5]
[mp13305-All-Series:10557] [11] simpleFoam() [0x41c5a0]
[mp13305-All-Series:10557] *** End of error message ***
--------------------------------------------------------------------------
mpirun noticed that process rank 0 with PID 10557 on node mp13305-All-Series exited on signal 8 (Floating point exception).

I guess it crashes because of the quality of the mesh? it is probably too big but with the inferior level of refinement the solution is not acceptable.

I have done your modifications for thefvSolution file.

And for sure the flow is turbulent (Re=100 000 in the channel), I use the realizable k-epsilon model.
Thanks for the help!

[Gerrit]

***Zero or negative cell volume detected.

This is as far as I know a no way. Could you figure out where this cells occur?
That makes it easier.

You can do that with the command

objToVTK -help

That should create a folder with the elements you quote in the command.

[MartinBlx]

Ok so the zeroVolumeCells check mesh should be the reason for the crash.

Here's there location. So it looks like a boundary issue right? Maybe the mesh is too fine?

https://drive.google.com/file/d/0B_T...ew?usp=sharing

I used the following code line to display the set in ParaView:

foamToVTK -cellSet zeroVolumeCells

EDIT: (I'm trying to mesh reducing the nCellsBetweenLevel feature in sHM. The checkMesh still fails but it seems to reduce the number of zeroVolumeCells, I only have 500 now. I guess I have to find the good value for not failing this check but still increase the resolution of my mesh compared to my first meshes. The CPU time is quite high but I have access to a far more powerful computer than before!)

[Gerrit]

that could solve the problem. By the way, how do you mesh your geometry? Do you use snappyHeymesh?

[MartinBlx]

SnappyHexMesh, implicit way (just one stl file with regions defined in it)