CFD Online Discussion Forums

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

psosnows October 6, 2011 16:10

Velocity/pressure/Temperature blow-up while dealing with unstructured meshes
3 Attachment(s)
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.

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

fvVectorMatrix UEqn
  + fvm::div(phi, U)
  - fvm::laplacian(mu/rho, U)
  + g * betaT * (T - T0)

The Temperature Equation:

volScalarField k = (K/(rho*Cp));
fvScalarMatrix TEqn
  + fvm::div(phi,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 consists of a 3D fluid box of dimension 0.5 x 0.5 x 0.5;
- front and back are cyclic.
- top, bottom, left and right patches have fixedValue (0 0 0) no-slip condition;
- 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;

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?).

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.

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.

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

kmooney October 18, 2011 09:14

Unstructured grids tend to be sensitive to your choice of discretization schemes.

Post #5 by Dr. Patterson in the thread below contains a few recommended schemes for skewed grids. You might want to give them a try if you haven't yet.

psosnows October 18, 2011 09:21

Thank you for the answer!
I am investigating this subject right now.

alberto October 19, 2011 14:48

Since it might interest others, there is a bug-report

alberto October 19, 2011 15:48

To put in practice what Henry suggested, use a limited scheme for the divergence. If your Pe is high, central differencing will give you problems.

Additionally, on unstructured grids, it is more appropriate to use leastSquares for gradients, since it preserves second-order accuracy.

You find a working set of schemes which fix the problem in the attachment.

All times are GMT -4. The time now is 17:32.