PIMPLE Stability

CallumG · July 11, 2022, 05:31

Morning Foamers,

I'm running into some stability issue when using interFOAM. My simulations crash at a similar time, consistently with a spike in "time step continuity errors" and the p_rgh residual. Here's the log of the last 2 timesteps computed:

Code:

Courant Number mean: 0.0230739546284 max: 38.197898741
Interface Courant Number mean: 0.000777623464901 max: 3.81051460639
Time = 14.88

PIMPLE: iteration 1
alpha.water BC on patch outlet
MULES: Solving for alpha.water
alpha.water BC on patch outlet
Phase-1 volume fraction = 0.83627896636  Min(alpha.water) = -103.911816209  Max(alpha.water) = 150.52296902
alpha.water BC on patch outlet
MULES: Solving for alpha.water
alpha.water BC on patch outlet
Phase-1 volume fraction = 0.836269635433  Min(alpha.water) = -487.723760213  Max(alpha.water) = 508.484315717
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.174349102868 )

Velocity BC on patch outlet
Force on turbine.blade1: 
(187539.627186 -54603.0964318 -9212.78760151)

Force on turbine.blade2: 
(196564.813807 20542.1016411 55167.9807565)

Force on turbine.blade3: 
(194943.762771 37863.6903348 -45806.9159198)

Force on turbine.hub: 
(5490.2140389 -17.5086780725 -2.9541132794)

Force on turbine.tower: 
(0 0 0)

Azimuthal angle (degrees) of turbine: 1250.42309209
Tip speed ratio of turbine: 4.4
Power coefficient from turbine: 0.400091257429
Rotor drag coefficient from turbine: 0.496560923074

DICPCG:  Solving for p_rgh, Initial residual = 0.0021608988582, Final residual = 0.000107723848319, No Iterations 137
time step continuity errors : sum local = 7.17961468072e-05, global = 2.95478051822e-06, cumulative = 4.04354720138e-05
DICPCG:  Solving for p_rgh, Initial residual = 0.000874326709373, Final residual = 4.33039666109e-05, No Iterations 55
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.174349102868 )

Velocity BC on patch outlet
time step continuity errors : sum local = 2.01227742579e-05, global = -2.21184730262e-06, cumulative = 3.82236247112e-05
DICPCG:  Solving for p_rgh, Initial residual = 0.000263194372774, Final residual = 9.86553797461e-08, No Iterations 246
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.174349102868 )

Velocity BC on patch outlet
time step continuity errors : sum local = 4.53912356615e-08, global = -1.17322367952e-09, cumulative = 3.82224514876e-05
smoothSolver:  Solving for omega, Initial residual = 0.974553445026, Final residual = 2.66883818386e-09, No Iterations 4
bounding omega, min: -116729.189533 max: 1019378.79645 average: 1.02686841445
smoothSolver:  Solving for k, Initial residual = 0.0188624168182, Final residual = 8.87519932106e-09, No Iterations 7
bounding k, min: -0.285459700159 max: 77.7821249136 average: 0.166107013953
ExecutionTime = 44122.8 s  ClockTime = 45602 s

Courant Number mean: 0.0299350494335 max: 2595.75221681
Interface Courant Number mean: 0.000886954295052 max: 492.131273028
Time = 14.885

PIMPLE: iteration 1
alpha.water BC on patch outlet
MULES: Solving for alpha.water
alpha.water BC on patch outlet
Phase-1 volume fraction = 0.836256700009  Min(alpha.water) = -25094.9966241  Max(alpha.water) = 58362.3120119
alpha.water BC on patch outlet
MULES: Solving for alpha.water
alpha.water BC on patch outlet
Phase-1 volume fraction = 0.836279047236  Min(alpha.water) = -9077126.4279  Max(alpha.water) = 14835489.7555
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.173966510369 )

Velocity BC on patch outlet
Force on turbine.blade1: 
(190766.74839 -56216.7106176 -9061.56576385)

Force on turbine.blade2: 
(199941.439777 21420.881685 56262.8625324)

Force on turbine.blade3: 
(196785.248964 38285.0426242 -47014.532404)

Force on turbine.hub: 
(5616.86786532 -28.3566064197 -4.57079845279)

Force on turbine.tower: 
(0 0 0)

Azimuthal angle (degrees) of turbine: 1250.84326114
Tip speed ratio of turbine: 4.4
Power coefficient from turbine: 0.410282425887
Rotor drag coefficient from turbine: 0.503842675799

DICPCG:  Solving for p_rgh, Initial residual = 0.128588696877, Final residual = 0.00624445596319, No Iterations 10
time step continuity errors : sum local = 0.00339774178095, global = -0.00075871394595, cumulative = -0.000720491494463
DICPCG:  Solving for p_rgh, Initial residual = 0.000102330247987, Final residual = 4.96600582123e-06, No Iterations 47
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.173966510369 )

Velocity BC on patch outlet
time step continuity errors : sum local = 0.00236246879967, global = -0.000305604529903, cumulative = -0.00102609602437
DICPCG:  Solving for p_rgh, Initial residual = 0.00172681281545, Final residual = 9.99392623372e-08, No Iterations 283
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.173966510369 )

Velocity BC on patch outlet
time step continuity errors : sum local = 3.61463698303e-06, global = 1.01674220237e-07, cumulative = -0.00102599435015
smoothSolver:  Solving for omega, Initial residual = 0.987600465769, Final residual = 5.33850547774e-10, No Iterations 5
bounding omega, min: -24776085.4187 max: 19488036.3905 average: 297.211121357
smoothSolver:  Solving for k, Initial residual = 0.323487262889, Final residual = 3.95184347366e-10, No Iterations 5
bounding k, min: -16.8027050907 max: 340.826867446 average: 0.173334322987
ExecutionTime = 44152.98 s  ClockTime = 45633 s

Courant Number mean: 0.305307922727 max: 132545.752862
Interface Courant Number mean: 0.0215072865315 max: 67758.6521143
Time = 14.89

PIMPLE: iteration 1
alpha.water BC on patch outlet
MULES: Solving for alpha.water
alpha.water BC on patch outlet
Phase-1 volume fraction = 0.836342658454  Min(alpha.water) = -121965323191  Max(alpha.water) = 109107175275
alpha.water BC on patch outlet
MULES: Solving for alpha.water
alpha.water BC on patch outlet
Phase-1 volume fraction = 0.836443830625  Min(alpha.water) = -2.88802783634e+15  Max(alpha.water) = 3.70574224922e+15
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.170576529572 )

Velocity BC on patch outlet
Force on turbine.blade1: 
(130484.732008 -30988.2423986 -4762.11120066)

Force on turbine.blade2: 
(145017.627717 14414.2680683 37040.5235279)

Force on turbine.blade3: 
(147517.277377 21830.9336265 -27213.8447653)

Force on turbine.hub: 
(3525.59840819 172.044707896 26.4389319004)

Force on turbine.tower: 
(0 0 0)

Azimuthal angle (degrees) of turbine: 1251.26343019
Tip speed ratio of turbine: 4.4
Power coefficient from turbine: 0.236732734963
Rotor drag coefficient from turbine: 0.362346920966

DICPCG:  Solving for p_rgh, Initial residual = 0.990013676426, Final residual = 0.0468017350473, No Iterations 9
time step continuity errors : sum local = 1488.54503318, global = -1.00783303658, cumulative = -1.00885903093
DICPCG:  Solving for p_rgh, Initial residual = 1.75086191517e-09, Final residual = 1.75086191517e-09, No Iterations 0
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.170576529572 )

Velocity BC on patch outlet
time step continuity errors : sum local = 46881.4838944, global = 3911.99295226, cumulative = 3910.98409323
DICPCG:  Solving for p_rgh, Initial residual = 6.42535826912e-07, Final residual = 9.55799719526e-08, No Iterations 3
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.170576529572 )

Velocity BC on patch outlet
time step continuity errors : sum local = 2559272.92531, global = -2526707.8432, cumulative = -2522796.85911
smoothSolver:  Solving for omega, Initial residual = 2.56912167053e-09, Final residual = 2.56912167053e-09, No Iterations 0
smoothSolver:  Solving for k, Initial residual = 0.703176894591, Final residual = 5.4883346337e-09, No Iterations 4
bounding k, min: -19.5825745916 max: 269.992183381 average: 0.17543350682
ExecutionTime = 44162.24 s  ClockTime = 45642 s

Courant Number mean: 44580578.5106 max: 3.53003684885e+13
Interface Courant Number mean: 968604.38205 max: 1.59633194071e+12
Time = 14.895

PIMPLE: iteration 1
alpha.water BC on patch outlet
MULES: Solving for alpha.water
alpha.water BC on patch outlet
Phase-1 volume fraction = -15937010.9633  Min(alpha.water) = -1.57713761489e+26  Max(alpha.water) = 2.9516660853e+26
alpha.water BC on patch outlet
MULES: Solving for alpha.water
alpha.water BC on patch outlet
Phase-1 volume fraction = -5.35687774577e+26  Min(alpha.water) = -3.63322160331e+50  Max(alpha.water) = 3.6404095442e+50
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.164308556659 )

Velocity BC on patch outlet
Force on turbine.blade1: 
(51631.2434405 -7756.81430588 -1133.86454448)

Force on turbine.blade2: 
(86866.2801199 7465.43665384 18775.4561849)

Force on turbine.blade3: 
(68488.8431049 4582.28635569 -5798.76745035)

Force on turbine.hub: 
(1779.66085343 254.583462175 37.214138432)

Force on turbine.tower: 
(0 0 0)

Azimuthal angle (degrees) of turbine: 1251.68359924
Tip speed ratio of turbine: 4.4
Power coefficient from turbine: 0.0797832917607
Rotor drag coefficient from turbine: 0.177345146484

DICPCG:  Solving for p_rgh, Initial residual = 0.0369707537481, Final residual = 0.00136598026531, No Iterations 8
time step continuity errors : sum local = 1.79053436895e+15, global = -12368703537, cumulative = -12371226333.8
DICPCG:  Solving for p_rgh, Initial residual = 4.68463524854e-24, Final residual = 4.68463524854e-24, No Iterations 0
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.164308556659 )

Velocity BC on patch outlet
time step continuity errors : sum local = 5.86444782243e+25, global = -2.76694252135e+25, cumulative = -2.76694252135e+25
DICPCG:  Solving for p_rgh, Initial residual = 0.999999978708, Final residual = 4.72823569472e-07, No Iterations 1000
3D_2D Absorption BC on patch inlet
"Correction Levels" 1( -0.164308556659 )

Velocity BC on patch outlet
time step continuity errors : sum local = 2.77996096845e+50, global = -7.60357064878e+48, cumulative = -7.60357064878e+48
[3] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[3] #1  Foam::sigFpe::sigHandler(int) at ??:?
[3] #2  ? in /lib64/libc.so.6
[3] #3  Foam::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, unsigned char, int) at ??:?
[3] #4  Foam::symGaussSeidelSmoother::scalarSmooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
[3] #5  Foam::symGaussSeidelSmoother::smooth(Foam::Field<double>&, Foam::Field<double> const&, unsigned char, int) const at ??:?
[3] #6  Foam::smoothSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[3] #7  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
[3] #8  Foam::fvMatrix<double>::solveSegregatedOrCoupled(Foam::dictionary const&) at ??:?
[3] #9  Foam::fvMesh::solve(Foam::fvMatrix<double>&, Foam::dictionary const&) const at ??:?
[3] #10  Foam::SolverPerformance<double> Foam::solve<double>(Foam::tmp<Foam::fvMatrix<double> > const&) at ??:?
[3] #11  Foam::kOmegaSSTStable<Foam::eddyViscosity<Foam::RASModel<Foam::IncompressibleTurbulenceModel<Foam::transportModel> > >, Foam::IncompressibleTurbulenceModel<Foam::transportModel> >::correct() at ??:?
[3] #12  ? at ??:?
[3] #13  __libc_start_main in /lib64/libc.so.6
[3] #14  ? at ??:?
[node1h30:26068] *** Process received signal ***
[node1h30:26068] Signal: Floating point exception (8)
[node1h30:26068] Signal code:  (-6)
[node1h30:26068] Failing at address: 0x162292000065d4
[node1h30:26068] [ 0] /lib64/libc.so.6(+0x36400)[0x2b25a4d31400]
[node1h30:26068] [ 1] /lib64/libc.so.6(gsignal+0x37)[0x2b25a4d31387]
[node1h30:26068] [ 2] /lib64/libc.so.6(+0x36400)[0x2b25a4d31400]
[node1h30:26068] [ 3] /exports/applications/apps/community/eng/OpenFOAM-1912/OpenFOAM-v1912/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam22symGaussSeidelSmoother6smoothERKNS_4wordERNS_5FieldIdEERKNS_9lduMatrixERKS5_RKNS_10FieldFieldIS4_dEERKNS_8UPtrListIKNS_17lduInterfaceFieldEEEhi+0x32d)[0x2b25a3bf99ed]
[node1h30:26068] [ 4] /exports/applications/apps/community/eng/OpenFOAM-1912/OpenFOAM-v1912/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam22symGaussSeidelSmoother12scalarSmoothERNS_5FieldIdEERKS2_hi+0x2d)[0x2b25a3bf9c2d]
[node1h30:26068] [ 5] /exports/applications/apps/community/eng/OpenFOAM-1912/OpenFOAM-v1912/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam22symGaussSeidelSmoother6smoothERNS_5FieldIdEERKS2_hi+0xd)[0x2b25a3bf968d]
[node1h30:26068] [ 6] /exports/applications/apps/community/eng/OpenFOAM-1912/OpenFOAM-v1912/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam12smoothSolver5solveERNS_5FieldIdEERKS2_h+0x5ac)[0x2b25a3bf12ac]
[node1h30:26068] [ 7] /exports/applications/apps/community/eng/OpenFOAM-1912/OpenFOAM-v1912/platforms/linux64GccDPInt32Opt/lib/libfiniteVolume.so(_ZN4Foam8fvMatrixIdE15solveSegregatedERKNS_10dictionaryE+0x11f)[0x2b259f79216f]
[node1h30:26068] [ 8] /exports/applications/apps/community/eng/OpenFOAM-1912/OpenFOAM-v1912/platforms/linux64GccDPInt32Opt/lib/libfiniteVolume.so(_ZN4Foam8fvMatrixIdE24solveSegregatedOrCoupledERKNS_10dictionaryE+0x27f)[0x2b259f0e6fef]
[node1h30:26068] [ 9] /exports/applications/apps/community/eng/OpenFOAM-1912/OpenFOAM-v1912/platforms/linux64GccDPInt32Opt/lib/libfiniteVolume.so(_ZNK4Foam6fvMesh5solveERNS_8fvMatrixIdEERKNS_10dictionaryE+0xf)[0x2b259f0a93bf]
[node1h30:26068] [10] /exports/applications/apps/community/eng/OpenFOAM-1912/OpenFOAM-v1912/platforms/linux64GccDPInt32Opt/lib/libincompressibleTurbulenceModels.so(_ZN4Foam5solveIdEENS_17SolverPerformanceIT_EERKNS_3tmpINS_8fvMatrixIS2_EEEE+0x126)[0x2b25a2e12ba6]
[node1h30:26068] [11] /home/s1453918/OpenFOAM/s1453918-v1912/platforms/linux64GccDPInt32Opt/lib/libturbulenceMultiphaseOlaFlowModels.so(_ZN4Foam15kOmegaSSTStableINS_13eddyViscosityINS_8RASModelINS_29IncompressibleTurbulenceModelINS_14transportModelEEEEEEES5_E7correctEv+0xd1a)[0x2b25bd66b71a]
[node1h30:26068] [12] olaFlow_current[0x444050]
[node1h30:26068] [13] /lib64/libc.so.6(__libc_start_main+0xf5)[0x2b25a4d1d555]
[node1h30:26068] [14] olaFlow_current[0x4490c5]
[node1h30:26068] *** End of error message ***
[node1h27][[33322,1],35][btl_tcp_frag.c:237:mca_btl_tcp_frag_recv] mca_btl_tcp_frag_recv: readv failed: Connection reset by peer (104)
[node1h27][[33322,1],34][btl_tcp_frag.c:237:mca_btl_tcp_frag_recv] mca_btl_tcp_frag_recv: readv failed: Connection reset by peer (104)
[node1h31][[33322,1],19][btl_tcp_frag.c:237:mca_btl_tcp_frag_recv] mca_btl_tcp_frag_recv: readv failed: Connection reset by peer (104)
[node1h27][[33322,1],32][btl_tcp_frag.c:237:mca_btl_tcp_frag_recv] mca_btl_tcp_frag_recv: readv failed: Connection reset by peer (104)
--------------------------------------------------------------------------
mpirun noticed that process rank 3 with PID 26068 on node node1h30 exited on signal 8 (Floating point exception).
--------------------------------------------------------------------------

and a copy of my fvSolution dictionary:

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "alpha.water.*"
    {
        nAlphaCorr      1;
        nAlphaSubCycles 2;
        alphaOuterCorrectors yes;
        cAlpha          1;

        MULESCorr       no;
        nLimiterIter    3;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0;
    }

    "pcorr.*"
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-5;
        relTol          0;
    }

    p_rgh
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-07;
        relTol          0.05;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-06;
        relTol          0;
    }

    "(k|epsilon|omega|B|nuTilda).*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-08;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor   no;
    nOuterCorrectors    1;
    nCorrectors         3;
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    fields
    {
    }
    equations
    {
        ".*" 1;
    }
}


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

It's also noteworthy that my timestep is fixed at 0.005 i.e. adjustTimeStep set to no. Also puzzling (to me) is that the line

Code:

Phase-1 volume fraction = XXXXXXX Min(alpha.water) = XXXXXXX  Max(alpha.water) = XXXXXXX

is consistently out of the [0,1] bounds, throughout the simulation, though when "stable" it's outside the bounds by a very tiny amount.

Does anyonne have any wisdom they could share on any of the above?

Cheers,
Callum

PassengerC07 · July 15, 2022, 05:20

You either have 1) a bad set of BCs, 2) bad quality mesh or 3) both 1) and 2).

Turn on the adjustTimeStep and see how the solver behaves (if the time step keeps getting smaller, think about whether the velocity in BC is reasonable). The mean courant number is not tiny compared to the max one so your mesh is probably not too fine throughout. Check the boundary layers and boundary conditions.

CallumG · July 15, 2022, 06:26

Hey Passenger,

I appreciate the reply. I expect the BC's are good because they are from olaFlow and run incredibly well in similar cases.

The mesh is essentially an empty rectangle that has refinement regions at the top and in the centre (I apply body forces there) with a non-orthogonality of <25, so I'm at pains to describe it as "bad quality" (I've attached a picture of a course version of the mesh).

An update on this is that I've successfully completed the simulation on a coarser mesh. Leading me to believe it's CFL related. So I'm currently in the process of testing at a lesser delta T (0.0025) which has surpassed the previous fail point. I am equally testing another version where I have upped maxCo to 2.0 and maxAlphaCo to 1.0. To hopefully make my simulations more efficient, I have enabled the adjustable time step for the new delta T, but I am limited as I need high(ish) frequency outputs.

All the best,
Callum

PassengerC07 · July 16, 2022, 04:17

Hi Callum,

The problem I have with CFL is similar to yours, that the BC and mesh quality is both good while the max Courant number is much larger than 1. The mean Courant number in my case is generally less than 0.01 but I still need to use adjustTimeStep to keep max Courant number at a low level. This gives me a time step at ~ 10^-5 s and it's a long simulation even for a case with 2s time span.

I think it's possible to pin point where the max Courant number occurs and adjust the mesh accordingly but I'm not there yet. My wild guess is the no-slip wall near the velocity inlet and the boundary layers are too fine.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
How to use PIMPLE properly?	floquation	OpenFOAM Running, Solving & CFD	26	June 18, 2024 05:38
PIMPLE – the value of the final under-relaxation factor	Zbynek	OpenFOAM	9	December 22, 2023 05:26
error while running modified pimple solver	R_21	OpenFOAM Programming & Development	0	May 28, 2015 06:59
A question on the PIMPLE algorithm	GerhardHolzinger	OpenFOAM Running, Solving & CFD	4	February 13, 2015 06:49
Stability in compressible pimple application	Thilo	OpenFOAM Programming & Development	0	May 30, 2012 11:30

July 15, 2022, 05:20		#2
PassengerC07 New Member Join Date: Jul 2022 Posts: 6 Rep Power: 4	You either have 1) a bad set of BCs, 2) bad quality mesh or 3) both 1) and 2). Turn on the adjustTimeStep and see how the solver behaves (if the time step keeps getting smaller, think about whether the velocity in BC is reasonable). The mean courant number is not tiny compared to the max one so your mesh is probably not too fine throughout. Check the boundary layers and boundary conditions.

July 16, 2022, 04:17		#4
PassengerC07 New Member Join Date: Jul 2022 Posts: 6 Rep Power: 4	Hi Callum, The problem I have with CFL is similar to yours, that the BC and mesh quality is both good while the max Courant number is much larger than 1. The mean Courant number in my case is generally less than 0.01 but I still need to use adjustTimeStep to keep max Courant number at a low level. This gives me a time step at ~ 10^-5 s and it's a long simulation even for a case with 2s time span. I think it's possible to pin point where the max Courant number occurs and adjust the mesh accordingly but I'm not there yet. My wild guess is the no-slip wall near the velocity inlet and the boundary layers are too fine.