Dear all,
I am running a simulation of a mixing tank on OF 7 where the only boundary condition that I have is Wall. No Inlet nor Outlet. I have a rotating zone for the impeller and stationary zone for the tank with the baffles. The geometry contains only half of the tank, due to the symmetry of the flow. The patches, where the tank is cut in half, communicate with each other via rotational cyclicAMI BC. The patches, where the rotating zone is connected to the stationary zone, communicate via a noOrdering cyclicAMI BC.
here is a checkMesh output:
Code:
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
Build : 7-109ba3c8d53a
Exec : checkMesh
Date : Aug 19 2019
Time : 13:46:04
Host : "fwd250"
PID : 12968
I/O : uncollated
Case : /mnt/sshfs/bigdata/openfoam/draw33/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/SSG/2007_Montante_onePhase_steadySSG
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time
Create polyMesh for time = 0
Time = 0
Mesh stats
points: 1134003
faces: 3298814
internal faces: 3197032
cells: 1082641
faces per cell: 6
boundary patches: 18
point zones: 0
face zones: 2
cell zones: 2
Overall number of cells of each type:
hexahedra: 1082641
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 0
Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
*Number of regions: 2
The mesh has multiple regions which are not connected by any face.
<<Writing region information to "0/cellToRegion"
<<Writing region 0 with 403920 cells to cellSet region0
<<Writing region 1 with 678721 cells to cellSet region1
Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
IMPELLER-SIDE_2 11880 12103 ok (non-closed singly connected)
BAFFLES 8712 9044 ok (non-closed singly connected)
TANK-WALL 11088 11571 ok (non-closed singly connected)
TANK-PERIODIC2 4620 4788 ok (non-closed singly connected)
TANK-PERIODIC1 4620 4788 ok (non-closed singly connected)
TANK-BOTTOM 3060 3216 ok (non-closed singly connected)
TANK-TOP 3060 3216 ok (non-closed singly connected)
IMPELLER-SIDE_1 11880 12103 ok (non-closed singly connected)
IMPELLER-BOTTOM 5108 5243 ok (non-closed singly connected)
IMPELLER-TOP 5108 5243 ok (non-closed singly connected)
IMPELLER-PERIODIC2 9435 9690 ok (non-closed singly connected)
IMPELLER-PERIODIC1 9435 9690 ok (non-closed singly connected)
IMPELLER3 1359 1392 ok (non-closed singly connected)
IMPELLER2 1359 1392 ok (non-closed singly connected)
IMPELLER1 1359 1392 ok (non-closed singly connected)
SHAFT 8000 8270 ok (non-closed singly connected)
SPARGER-SURFACE 1200 1274 ok (non-closed singly connected)
INLET 499 548 ok (non-closed singly connected)
Checking geometry...
Overall domain bounding box (-0.118 -0.118 0) (0.118 1.44504e-17 0.236)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (1.3176e-16 1.39255e-15 -4.29463e-17) OK.
Max cell openness = 4.42129e-16 OK.
Max aspect ratio = 29.1252 OK.
Minimum face area = 6.71041e-08. Maximum face area = 7.10408e-05. Face area magnitudes OK.
Min volume = 3.7373e-11. Max volume = 1.51012e-07. Total volume = 0.00510649. Cell volumes OK.
Mesh non-orthogonality Max: 45.2912 average: 10.76
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.920013 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
The impeller is rotating with 450 rpm.
First I ran the simulation with kOmegaSST. I did not get the best convergence for the pressure but at least the simulation was stable. When I am running the SSG model, the simulation first runs normally then suddenly diverges and crashes. I noticed a
bounding epsilon message which then leads to a very big
time step continuity error. This leads afterwards to very small residual in velocity (1e-19 e.g.) and big one in pressure (1 e.g.). These residuals then oscillates then the whole simulation stops.
I started then reading about these kind of issues in the forum. All changes of numerical schemes in fvSchemes or the solvers in fvSolution did not really help too much.
fvSchemes:
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | cfMesh: A library for mesh generation |
| \\ / O peration | |
| \\ / A nd | Author: Franjo Juretic |
| \\/ M anipulation | E-mail: franjo.juretic@c-fields.com |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default steadyState;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,U) bounded Gauss limitedLinearV 1;
div(phi,k) Gauss upwind; //bounded Gauss limitedLinear 1;
div(phi,epsilon) Gauss upwind; //bounded Gauss limitedLinear 1;
div(phi,omega) bounded Gauss limitedLinear 1;
div(phi,R) bounded Gauss limitedLinear 1; //upwind;
div((nuEff*dev2(T(grad(U))))) Gauss linear;
div((nu*dev2(T(grad(U))))) Gauss linear;
div(R) Gauss linear;
}
laplacianSchemes
{
default Gauss linear uncorrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default uncorrected;
}
fluxRequired
{
default no;
p;
}
wallDist
{
method Poisson;
nRequired true;
}
// ************************************************************************* //
fvSolution:
Code:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | cfMesh: A library for mesh generation |
| \\ / O peration | |
| \\ / A nd | Author: Franjo Juretic |
| \\/ M anipulation | E-mail: franjo.juretic@c-fields.com |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
"(k|epsilon|omega|R).*"
{
solver PBiCG;
//solver smoothSolver;
preconditioner DILU;
tolerance 1e-8;
relTol 0.001;
minIter 1;
}
yPsi
{
solver PCG;
preconditioner none;
tolerance 1e-10;
relTol 0;
}
p
{
solver GAMG;
tolerance 1e-8;
relTol 1e-8;
minIter 5;
maxIter 100;
smoother GaussSeidel;
nPreSweeps 1;
nPostSweeps 3;
nFinestSweeps 3;
scaleCorrection true;
directSolveCoarsest false;
cacheAgglomeration on;
nCellsInCoarsestLevel 1000;
agglomerator faceAreaPair;
mergeLevels 1;
minIter 2;
};
U
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-09;
relTol 0;
};
}
SIMPLE
{
nNonOrthogonalCorrectors 4;
pRefCell 0;
pRefValue 0;
residualControl
{
p 1e-5;
U 1e-5;
"(k|epsilon|omega|R)" 1e-5;
}
}
relaxationFactors
{
fields
{
p 0.3;
}
equations
{
U 0.5;
k 0.3;
epsilon 0.3;
omega 0.3;
R 0.3;
}
}
PISO
{
nCorrectors 5;
nNonOrthogonalCorrectors 2;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //
I tried running kEpsilon model first then mapping the results to the SSG model but this also did not help. I started plotting the residuals of all the variables including epsilon and reynolds stresses, and then I saw that after a while of running the simulation,
the value of the residual of both epsilon and reynolds stresses start to increase. Then the same happens to the residual of velocity, while the residual of pressure is decreasing in value untill the whole simulation crashes.
So I started playing around with the relaxation factors. When setting the under relaxation factor of velocity and epsilon equal to 0.001 and reynolds stresses to 0.0001, the simulation is much more stable. The residual of reynolds stresses still increase in value but much more slower than before. After running the simulation with these relaxation factors for a while, the pressure residual starts to oscillates around 0.0002 and continues like this.
Code:
relaxationFactors
{
fields
{
p 0.3;
}
equations
{
U 0.001;
k 0.3;
epsilon 0.001;
omega 0.3;
R 0.0001;
}
}
This sounds to me like I am stopping everything from being solved expect for the pressure, right? If I reach a "good" convergence with this set-up by either monitoring the residual or convergence probes, would the solution be of SSG model? Is there anything else I should pay attention to? Did I do anything wrong setting up this case? Why are RSTM much more sensetive than 2-equation models?
I want to map the best converged solution of the steady-state case to the transient case and start again. When I run the transient case now, I am getting oscillations faster and the simulation crashes also quicker than in steady-state case.