[sebbi]

Hi,

I'm trying to simulate a dust enclosure. The first I want to simulate is the closed enclosure. I have one outlet (with a given velocity of about 20 m/s) and 2 or 3 inlets (with a given pressure of 0 Pa). Later I want to simulate with open door (pressure 0 Pa), to look if there is an outflow. For this I had made a model in *.stl. I meshed this model using snappyHexMesh. Now I have the problem, that the solution does not converge. Sometimes it blows up, and sometimes the residuals stay at one point and go up and down slowly (many iterations). So I can't get a right solution. I also viewed the results in paraview, but I can't see a steady solution. For example:

Time = 815

DILUPBiCG:  Solving for Ux, Initial residual = 0.00381937156, Final residual = 5.46407448e-07, No Iterations 2
DILUPBiCG:  Solving for Uy, Initial residual = 0.0039635767, Final residual = 3.38204802e-07, No Iterations 2
DILUPBiCG:  Solving for Uz, Initial residual = 0.00394110401, Final residual = 5.52425189e-07, No Iterations 2
DICPCG:  Solving for p, Initial residual = 0.0339215183, Final residual = 3.39079191e-06, No Iterations 260
time step continuity errors : sum local = 1.59943344e-05, global = 6.05402128e-10, cumulative = 1.90916487e-06
ExecutionTime = 3293.45 s  ClockTime = 3304 s (has to be about 10 000sec I made a break)

My checkMesh now looks like this:

Mesh stats
    points:           590512
    faces:            1672734
    internal faces:   1609959
    cells:            541732
    boundary patches: 8
    point zones:      0
    face zones:       0
    cell zones:       0

Overall number of cells of each type:
    hexahedra:     516514
    prisms:        10178
    wedges:        145
    pyramids:      0
    tet wedges:    322
    tetrahedra:    0
    polyhedra:     14573

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
    defaultFaces        0        0        ok (empty)
    xxxxxx_lalu          2179     2361     ok (non-closed singly connected)
    xxxxxx_outlet        271      311      ok (non-closed singly connected)
    xxxxxx_spa_links  2271     2387     ok (non-closed singly connected)
    xxxxxx_spa_rechts2270     2386     ok (non-closed singly connected)
    xxxxxx_tuerschl      131      249      ok (non-closed singly connected)
    xxxxxx_tuer          3168     3337     ok (non-closed singly connected)
    xxxxxx_walls         52485    54205    ok (non-closed singly connected)

Checking geometry...
    Overall domain bounding box (-0.0495049591 -1.879 0) (1.1265 0 1.87)
    Mesh (non-empty, non-wedge) directions (1 1 1)
    Mesh (non-empty) directions (1 1 1)
    Boundary openness (-2.18009372e-18 3.97104586e-19 2.29753118e-19) OK.
    Max cell openness = 2.59771373e-16 OK.
    Max aspect ratio = 37.4175372 OK.
    Minumum face area = 1.79977657e-06. Maximum face area = 0.00207470849.  Face area magnitudes OK.
    Min volume = 2.31921363e-09. Max volume = 7.75293658e-05.  Total volume = 2.47929553.  Cell volumes OK.
    Mesh non-orthogonality Max: 64.3808597 average: 5.39379502
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 3.82309381 OK.

Mesh OK.

For me it looks good. So I made a set up for a laminar incompressible flow. (I think, the effect of compression is very little and laminar is a good choice for beginnig). My parameters are:

RASModell laminar
turbulence off
simulationType laminar
nu=13.5e-06 (constant for air)

For the solver I choose the following:

p
{
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-07;
        relTol          0.001;}
U
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-05;
        relTol          0.001;
    }
SIMPLE
{
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    p               0.3;
    U               0.2;
}

I also tried lower and higher relaxationFactors GAMG for p and U, but nothing helps. I also use the following schemes:

div(phi,U)      Gauss upwind;
interpolate(U)  linear;
ddt default steadyState;

The BC is like this:

U
outlet -> fixedValue uniform (0 0 20)
inlets -> zeroGradient
walls and doors -> fixedValue uniform (0 0 0)

p
outlet -> zeroGradient
inlets -> fixedValue 0
walls and doors -> zeroGradient

Now I tried to get rid of the problems for two weeks. I tried other BCs (outletInlet), other schemes (linear instead of upwind) and so on.

I hope you can help me with this problem and give me an hint to clear the fault I surely made somewhere. If I should give you further informations don't hesitate to ask.

Thank you very much.

Sebastian

[linnemann]

Hi

your skewness is very bad (3.8), best results coming from having it below 0.9. you could also try leastSquares (or faceLimited) as gradSchemes and maybe gammaV 0.2 for div(phi,U). and upwind for the turbulence properties. I would still focus on the grid skewness as my main place of interest. Try using another mesh generator like Netgen, Engrid or Salome.

ddtSchemes
{
    default steadyState;
}

gradSchemes
{
    default         leastSquares;
    grad(p)         leastSquares;
    grad(U)         leastSquares;
//    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss gammaV 0.2;
    div(phi,k)      Gauss upwind;
    div(phi,omega)  Gauss upwind;
    div((nuEff*dev(grad(U).T()))) Gauss linear;
}

laplacianSchemes
{
 //   default         Gauss linear corrected;
    default         Gauss linear limited 0.5;
//    default         Gauss linear limited 0.333;
}

interpolationSchemes
{
    default         linear;
    interpolate(U)  linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    p;
}

// ************************************************************************* //

Also using smoothSolver for U,k,omega,epsilon would be an idea same as using GAMG for p.

 
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;
    };

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

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

SIMPLE
{
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    p               0.3;
    U               0.7;
    k               0.7;
    omega           0.7;
}

// ************************************************************************* //

[linnemann]

Hi again

another thing if you have it write out values often you can open paraview and examine where the solution values blow up. This way you can examine the mesh where this occurs and I bet it is where you have your most skewed cells.

I once had one badly skewed cell (4.12) and it made the whole solution crash. So even though checkMesh states the mesh is ok you have to use your common sense to see if "you" accept it to be ok :-).

[sebbi]

Yes you're right, I saw this phenomenon before, that two skewed cells at the inlet make trouble. I refinde the mesh using snappyHexMesh with refinementRegion or sth. like that. The mesh I had bevor was even worse than this.
I decided to use sHM cause, the other tools often crash (netgen) or do something for long time without feedback (Salome). But I think I have no other chance and so I'll give Salome a second chance and try it.
The schemes you told me, I will try them tomorrow (probably on the new mesh).

Thank you for your help