 |
 |
 |
|
|
|
[Sponsors] | ||||
March 5, 2024, 17:29
|
Error when using Crank Nicolson ddt scheme
|
#1 |
|
Member
Pedro Gouveia
Join Date: Oct 2022
Location: Portugal
Posts: 59
Rep Power: 3  |
Hey guys,
I am performing a centrifugal pump simulation in OF11. At the moment I am trying a transient approach, since I have already ran all the steady state simulations I needed. The transient simulation is started with the results from the steady state simulation, since I do not care for the startup of the pump. It runs at 620 rpm, with 0.0063m^3/s flow rate. It is a relatively small pump. The simulation is dynamic, using the sliding mesh approach and the new Non Conformal Coupling feature from the foundation branch. I am using vanLeer schemes for turbulence (using kEpsilon model) and LUST for convective term. No cell/face limiters are being used (although I already tried and it didn't improve anything). When using Euler for the ddt schemes, no error appears. However, when using Crank Nicolson 0.1 to 0.9 a floating point exception occurs. Code:
-------------------------------------------------------------------------- Primary job terminated normally, but 1 process returned a non-zero exit code. Per user-direction, the job has been aborted. -------------------------------------------------------------------------- -------------------------------------------------------------------------- mpirun noticed that process rank 3 with PID 0 on node omnideayoga-Yoga-Slim-7-Pro-14ARH7 exited on signal 8 (Floating point exception). Code:
Mesh stats
points: 801625
faces: 2283864
internal faces: 2148432
cells: 741672
faces per cell: 5.97609
boundary patches: 30
point zones: 0
face zones: 3
cell zones: 1
Overall number of cells of each type:
hexahedra: 723936
prisms: 17736
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: 3
The mesh has multiple regions which are not connected by any face.
<<Writing region information to "0/cellToRegion"
<<Writing region 0 with 335968 cells to cellSet region0
<<Writing region 1 with 246432 cells to cellSet region1
<<Writing region 2 with 159272 cells to cellSet region2
Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
NCC3_2 6976 7412 ok (non-closed singly connected)
wall_3_2 1680 1802 ok (non-closed singly connected)
wall_3_3 2000 2142 ok (non-closed singly connected)
wall_3_4 2240 2397 ok (non-closed singly connected)
wall_3_layers 28276 28964 ok (non-closed singly connected)
wall_tongue_1 128 153 ok (non-closed singly connected)
wall_3_1 1360 1462 ok (non-closed singly connected)
wall_tongue_2 128 153 ok (non-closed singly connected)
wall_tongue_3 160 187 ok (non-closed singly connected)
outlet 784 850 ok (non-closed singly connected)
wall_4 4480 4794 ok (non-closed singly connected)
wall_4_layers 13720 14100 ok (non-closed singly connected)
NCC2_1 2432 2584 ok (non-closed singly connected)
NCC2_3 6976 7412 ok (non-closed singly connected)
blade_LE 320 408 ok (non-closed singly connected)
blade_pressure 7808 8364 ok (non-closed singly connected)
blade_suction 7936 8500 ok (non-closed singly connected)
blade_tip 448 544 ok (non-closed singly connected)
wall_2 30804 31908 ok (non-closed singly connected)
inlet 1852 1899 ok (non-closed singly connected)
NCC1_2 2432 2584 ok (non-closed singly connected)
wall_1 12492 12691 ok (non-closed singly connected)
nonConformalCouple_12_on_NCC1_20 0 ok (empty)
nonConformalCouple_12_on_NCC2_10 0 ok (empty)
nonConformalCouple_23_on_NCC2_30 0 ok (empty)
nonConformalCouple_23_on_NCC3_20 0 ok (empty)
nonConformalError_on_NCC2_30 0 ok (empty)
nonConformalError_on_NCC1_20 0 ok (empty)
nonConformalError_on_NCC3_20 0 ok (empty)
nonConformalError_on_NCC2_10 0 ok (empty)
Checking geometry...
Overall domain bounding box (-0.153983 -0.18447 -0.0123) (0.220992 0.482614 0.1647)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (-3.87232e-17 -5.16145e-18 -6.18043e-16) OK.
Max cell openness = 3.32906e-16 OK.
Max aspect ratio = 13.8188 OK.
Minimum face area = 9.82531e-08. Maximum face area = 1.4498e-05. Face area magnitudes OK.
Min volume = 1.23456e-10. Max volume = 2.82255e-08. Total volume = 0.00389344. Cell volumes OK.
Mesh non-orthogonality Max: 54.0036 average: 7.95413
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 2.97155 OK.
Coupled point location match (average 0) OK.
Mesh OK.
End
Code:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 11
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default CrankNicolson 0.5;
}
gradSchemes
{
default Gauss linear;
}
divSchemes
{
default none;
div(phi,k) bounded Gauss vanLeer;
div(phi,epsilon) bounded Gauss vanLeer;
div(phi,U) bounded Gauss LUST grad(U);
div((nuEff*dev2(T(grad(U))))) bounded Gauss linear;
}
laplacianSchemes
{
default Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
wallDist
{
method Poisson;
}
// ************************************************************************* //
Code:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 11
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver GAMG;
tolerance 1e-08;
relTol 0.01;
smoother GaussSeidel;
nCellsInCoarsestLevel 1000;
//nSweeps 200;
minIter 5;
}
pFinal
{
$p;
tolerance 1e-08;
relTol 0.0;
}
"pcorr.*"
{
$p;
tolerance 1e-03;
relTol 0.0;
}
U
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-08;
//nSweeps 200;
relTol 0.01;
//type coupled;
//solver PBiCGStab;
//preconditioner DILU;
//tolerance (1e-08 1e-08 1e-08);
//tolerance 1e-14
//relTol (0 0 0);
//relTol 0.0;
minIter 3;
}
UFinal
{
$U;
tolerance 1e-10;
relTol 0.0;
}
"(k|epsilon|omega)"
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-08;
relTol 0.01;
//nSweeps 50;
minIter 2;
}
"(k|epsilon|omega)Final"
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-10;
relTol 0.0;
//nSweeps 50;
minIter 2;
}
MeshPhi
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-2;
relTol 0.0;
}
}
potentialFlow
{
nNonOrthogonalCorrectors 10;
}
PIMPLE
{
correctMeshPhi yes;
momentumPredictor yes;
nNonOrthogonalCorrectors 2;
nCorrectors 2;
nOuterCorrectors 30;
residualControl
{
p 1e-8;
U 1e-10;
k 1e-8;
epsilon 1e-8;
}
}
relaxationFactors
{
fields
{
p 0.3;
pFinal 1;
}
equations
{
"(U|k|omega)" 0.7;
"(U|k|omega)Final" 1;
}
}
// ************************************************************************* //
|
|
|
|
|
|
March 6, 2024, 07:22
|
|
#2 |
|
Senior Member
Jan
Join Date: Jul 2009
Location: Hamburg
Posts: 137
Rep Power: 19 |
Hi Pedro,
I have just checked with a test case based on the propeller tutorial and can reproduce the behaviour. For a case with moving mesh (solid body rotation) in combination with NCC the sovler crashes at second time step as soon as one uses CrankNicolson ddtScheme. I think, this is most likely due to the NCC, and there are some values (fluxes or whatever) required from the previous time step, but are not available. Looks like this is eihter a bug or the combination of NCC and secondOrder time integration is not (yet) implemented. If you really require CrankNicolson, you may prepare a simple test case and submit a bug report. Btw, any reason for using such huge amount of outer correctors? Best, Jan |
|
|
|
|
|
|
March 6, 2024, 08:02
|
|
#3 | |
|
Member
Pedro Gouveia
Join Date: Oct 2022
Location: Portugal
Posts: 59
Rep Power: 3 |
Quote:
Do you know how can I submit a bug report? I tried with the backward scheme and everything works fine. The reason I am using this many outer correctors is because I read that we should not limit the pimple iterations, but use a convergence criteria so that the cycle exits automatically before the max number of outer iterations. However, when I do that, nothing really happens. I am using this fvSolution: Code:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 11
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver GAMG;
tolerance 1e-06;
relTol 0.1;
smoother GaussSeidel;
nCellsInCoarsestLevel 1000;
}
pFinal
{
$p;
relTol 0.0;
}
"pcorr.*"
{
$p;
tolerance 1e-03;
relTol 0.0;
}
U
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-06;
relTol 0.1;
}
UFinal
{
$U;
relTol 0.0;
}
"(k|epsilon|omega)"
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-06;
relTol 0.01;
}
"(k|epsilon|omega)Final"
{
solver smoothSolver;
smoother GaussSeidel;
relTol 0.0;
//nSweeps 50;
minIter 2;
}
MeshPhi
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-2;
relTol 0.0;
}
}
potentialFlow
{
nNonOrthogonalCorrectors 10;
}
PIMPLE
{
correctMeshPhi yes;
momentumPredictor yes;
nNonOrthogonalCorrectors 2;
nCorrectors 2;
nOuterCorrectors 30;
residualControl
{
p 1e-5;
U 1e-6;
//k 1e-8;
//epsilon 1e-8;
//".*" 1e-5;
}
}
relaxationFactors
{
fields
{
p 0.4;
pFinal 1;
}
equations
{
"(U|k|omega)" 0.4;
"(U|k|omega)Final" 1;
}
}
// ************************************************************************* //
|
||
|
|
|
||
|
March 6, 2024, 08:29
|
|
#4 |
|
Senior Member
Jan
Join Date: Jul 2009
Location: Hamburg
Posts: 137
Rep Power: 19 |
Ok, interessting. If it runs with backward scheme, it seems to be a specific problem of CrankNicolson, not of 2nd order time schemes in general.
As the backward scheme is much more stable (with similar accuracy than CN) , it suggest to use this is your simulations. Bugs can be reported via: https://bugs.openfoam.org/rules.php I would not use the residual control, and just set the outer correctors to a fixed number (3 to 5 is usually fine, but depend on the case (mesh, time step, flow type) and the piso ( nCorrectors) to 1. |
|
|
|
|
|
|
March 6, 2024, 08:31
|
|
#5 | |
|
Member
Pedro Gouveia
Join Date: Oct 2022
Location: Portugal
Posts: 59
Rep Power: 3 |
Quote:
|
||
|
|
|
||
 |
|
|
Similar Threads
|
||||
| Thread | Thread Starter | Forum | Replies | Last Post |
| Second-order backward differencing vs Crank Nicolson | Beans8 | Main CFD Forum | 7 | June 23, 2023 04:59 |
| Implementtaion of Crank Nicolson scheme in OpenFOAM | hchen | OpenFOAM Running, Solving & CFD | 5 | May 6, 2020 14:25 |
| Crank Nicolson scheme implemented wrong? | rajibroy | OpenFOAM Programming & Development | 10 | May 5, 2020 09:57 |
| Crank Nicolson scheme | msrinath80 | OpenFOAM Running, Solving & CFD | 6 | November 14, 2010 13:59 |
| Crank Nicolson scheme contd | msrinath80 | OpenFOAM Running, Solving & CFD | 0 | December 12, 2006 22:41 |
Linear Mode
