CFD Online Discussion Forums - Velocity/pressure/Temperature blow-up while dealing with unstructured meshes

Good day!

We would like to present here a very disturbing issue regarding solving cases with unstructured meshes with OpenFOAM.
We are holding from calling it a "major bug" just because we need someone to confirm the presented results.
The issue arises with usage of non-regular meshes for solving buoyant problems.

At certain moment during a simulation, the pressure, temperature and velocity field blows up. This always starts in a single, random, unstructured element within the fluid.

THE SOLVER
For presentation of the bug, we use an incompressible PISO-based, transient, buoyant solver (attached to this thread in the .zip).
The Velocity equation:

Code:

fvVectorMatrix UEqn

(

    fvm::ddt(U)

  + fvm::div(phi, U)

  - fvm::laplacian(mu/rho, U)

  + g * betaT * (T - T0)

);

The Temperature Equation:

Code:

volScalarField k = (K/(rho*Cp));

k.correctBoundaryConditions();

fvScalarMatrix TEqn

(

    fvm::ddt(T)

  + fvm::div(phi,T)

 ==

    fvm::laplacian(k,T)

);

where U is a volVectorField, phi is standard-calculated, mu, rho, Cp, T, K are volScalarFields, g is a dimmensionedVector and T0 is a dimmensionedScalar.

THE TEST CASE
Geometry:
- the test case consists of a 3D fluid box of dimension 0.5 x 0.5 x 0.5;
- front and back are cyclic.
Velocity:
- top, bottom, left and right patches have fixedValue (0 0 0) no-slip condition;
Temperature:
- top and bottom are zeroGradient;
- left has fixed value 1;
- right has fixed value 0;
- for each new mesh, the steady, linear distribution of temperature field was obtained using laplacianFoam for initial condition;
Other fields' amd properties set-up:
- top, bottom, left and right are zeroGradient;
- density rho [1 -3 0 0 0 0 0] 1000;
- specific heat Cp [0 2 -2 -1 0 0 0] 1000;
- dynamic visc. mu [1 -1 -1 0 0 0 0] 0.1;
- thermal cond. K [1 1 -3 -1 0 0 0] 1;
- therm. expansion coeff betaT [0 0 0 -1 0 0 0] 0.01;
- reference temp. T0 [0 0 0 1 0 0 0] 0.5;
- gravity acc. g [0 1 -2 0 0 0 0] ( 0 -9.81 0 );
- max Courant number maxCoNum 0.3;
Schemes and precisions can be found in the attached cases.

All test cases are attached in the .zip

Steps to reproduce:
- unpack ;)
- blockMesh;
- for "Pyramid" case, run modifyMesh;
- change fvSolution to allow running laplacianFoam;
- run laplacianFoam to acquire linear temperature distribution;
- change back fvSolution;
- use obtained linear T field as initial condition for attached controlBuoyantFoam solver;
- run the cases;

REFERENCE RUN
The reference case was meshed using blockMesh, with regular mesh with number of cells 40 x 40 x 40. The mesh is not fine enough near the walls and will produce errors near the walls, but considering problems presented later, those errors can be omitted (but a question arises- why do they occur anyway?).

CROSS 2x2x2 REFINEMENT
The mesh was refined in single cell lines in 3D, to acquire splitting of one cell into 2x2x2 hexes (see hex2x2x2Mesh.jpg). The simulation was performed giving same results as Reference run.

PIRAMID SPLITTING
One cell of Reference mesh was was split into 6 identical, regular pyramids using modifyMesh tool. The simulation was performed, and after working for a while, the simulation broke. Cells near the refinement developed non-physical large temperature and velocity (see TempThreshold.jpg).
After initial test, the mesh was treated with renumberMesh tool, but the simulation gave same result.

CONCLUSIONS AND REMARKS
The presented cases were designed to visualise the error. We encountered this strange behaviour many times, always while using unstructured mesh. The simulations blew up in random mesh positions for no particular reasons. We have run them with increased number of PISO and non-orthogonal corrections (even up to 16 PISO and 8 non-ort). And we always encountered problems.

Industrial applications in most cases require unstructured grids. It is true, that maybe one can play with the mesh and pray that this time the simulation will not blow up... but if a simple solver does not handle simple mesh, that means something is wrong.

Right now we are investigating the "deep" OF code checking for the reason of the problem. We are looking forward to start some collaboration to eliminate it.

Anxious for Your kind replays,
Pawel S and Colleagues