[TurbulentGuy]

Hello all
I am trying to simulate a cavitating flow inside a venturi-section using InterPhaseChangeFoam- an unsteady, multiphase solver.
Here is my fvSolution

/*--------------------------------*- 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;
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

        MULESCorr       yes;
        nLimiterIter    5;

        solver          PBiCG;
	preconditioner	DILU;
        tolerance       1e-8;//1e-8
        relTol          0;
    };

    "U.*"
    {
        solver           PBiCG;
        preconditioner	 DILU;
        tolerance        1e-6; //1e-6
        relTol           0;
    };

    "(k|omega)"
    {
        solver           PBiCG;
        preconditioner	 DILU;
        tolerance        1e-6; //1e-6
        relTol           0;
    };
 

    "(U|k|omega)Final"
    {
        solver           PBiCG;
	preconditioner	 DILU;
	tolerance	 1e-6; //1e-6
        relTol 		 0;
    }
    p_rgh
    {
        solver           GAMG;
        tolerance        1e-8; //1e-8
        relTol           0.1;

        smoother         DICGaussSeidel;
        nPreSweeps       0;
        nPostSweeps      2;

        cacheAgglomeration true;

        nCellsInCoarsestLevel 10;
        agglomerator     faceAreaPair;
        mergeLevels      1;

        maxIter          50;
    };

    nuTilda
    {
	solver		PBiCG;
	preconditioner	DILU;
	smoother	DICGaussSeidel;
	nSweeps		2;
	tolerance	1e-6;
	relTol		0;
    };

    p_rghFinal
    {
        solver            PCG;
	preconditioner
        {
            preconditioner  GAMG;

            tolerance        1e-6; //1e-6
            relTol           0;

            nVcycles         2;

            smoother         DICGaussSeidel;
            nPreSweeps       0;
            nPostSweeps      2;
            nFinestSweeps    2;

            cacheAgglomeration true;
            nCellsInCoarsestLevel 10;
            agglomerator     faceAreaPair;
            mergeLevels      1;
        };
        tolerance        1e-7; //1e-7
        relTol           0;
        maxIter          50;
    };
}

PIMPLE
{
    momentumPredictor          yes;
    nOuterCorrectors           100;
    nCorrectors                3;
    nNonOrthogonalCorrectors   0;
    residualControl  //outerCorrectorResidualControl
    	{
        	"(p_rgh|U)"
		{
		 tolerance	 5e-5;//5e-5
		 relTol 	 0;//0.1
		}
		"(k|omega|epsilon)"
		{
		 tolerance	1e-4;
		 relTol		0;		
		}
    	}	
}


relaxationFactors
{
    fields
    {
	p		0.8;
    }
    equations
    {
        U               0.8; // 0.9 is more stable but 0.95 more convergent
        ".*"            0.8; // 0.9 is more stable but 0.95 more convergent
	k		0.8;
	omega		0.8;
    }
}

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

As you can see I have applied a residual control limit of 5e-5 on pressure and velocity and 1e-4 on k and omega(I am using k-omega SST model here).
However, looking at the residual files, the iterations continue even after dropping to the required level

Courant Number mean: 0.0528581 max: 2.48597
Time = 1.39901

PIMPLE: iteration 1
DILUPBiCG:  Solving for alpha.water, Initial residual = 3.68011e-07, Final residual = 1.36401e-09, No Iterations 1
Phase-1 volume fraction = 0.991636  Min(alpha.water) = 0.0154223  Max(alpha.water) = 1
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Liquid phase volume fraction = 0.991616  Min(alpha.water) = 0.0154223  Max(alpha.water) = 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.00127737, Final residual = 7.84762e-07, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.00227295, Final residual = 3.61585e-09, No Iterations 2
GAMG:  Solving for p_rgh, Initial residual = 6.75951e-06, Final residual = 5.29966e-10, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 1.63286e-07, Final residual = 3.04435e-10, No Iterations 1
GAMGPCG:  Solving for p_rgh, Initial residual = 1.90892e-08, Final residual = 1.90892e-08, No Iterations 0
This is kOmegaSSTF
DILUPBiCG:  Solving for omega, Initial residual = 0.000120964, Final residual = 7.63802e-08, No Iterations 1
bounding omega, min: -29.5739 max: 2.83886e+06 average: 82685.8
DILUPBiCG:  Solving for k, Initial residual = 0.00256498, Final residual = 2.6164e-09, No Iterations 2
bounding k, min: -0.100634 max: 22.1979 average: 0.703296
PIMPLE: iteration 2
DILUPBiCG:  Solving for alpha.water, Initial residual = 1.51834e-07, Final residual = 3.8706e-11, No Iterations 1
Phase-1 volume fraction = 0.991644  Min(alpha.water) = 0.0154281  Max(alpha.water) = 1
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Liquid phase volume fraction = 0.991624  Min(alpha.water) = 0.0154281  Max(alpha.water) = 1
DILUPBiCG:  Solving for Ux, Initial residual = 0.000288568, Final residual = 1.69436e-07, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.000502123, Final residual = 2.93426e-07, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 6.27618e-07, Final residual = 6.02198e-10, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 2.07927e-07, Final residual = 3.39953e-10, No Iterations 1
GAMGPCG:  Solving for p_rgh, Initial residual = 1.89539e-08, Final residual = 1.89539e-08, No Iterations 0
This is kOmegaSSTF
DILUPBiCG:  Solving for omega, Initial residual = 3.04491e-05, Final residual = 1.35276e-08, No Iterations 1
bounding omega, min: -22.4256 max: 2.83873e+06 average: 82690.8
DILUPBiCG:  Solving for k, Initial residual = 0.000716519, Final residual = 5.31525e-07, No Iterations 1
bounding k, min: -0.0272711 max: 22.274 average: 0.703536
PIMPLE: iteration 3
DILUPBiCG:  Solving for alpha.water, Initial residual = 1.46881e-07, Final residual = 3.71256e-11, No Iterations 1
Phase-1 volume fraction = 0.991645  Min(alpha.water) = 0.0154305  Max(alpha.water) = 1
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Liquid phase volume fraction = 0.991625  Min(alpha.water) = 0.0154305  Max(alpha.water) = 1
DILUPBiCG:  Solving for Ux, Initial residual = 6.57397e-05, Final residual = 3.9333e-08, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 0.000114833, Final residual = 7.40035e-08, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 1.85128e-07, Final residual = 3.77726e-10, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 2.61274e-08, Final residual = 2.8472e-10, No Iterations 1
GAMGPCG:  Solving for p_rgh, Initial residual = 1.89681e-08, Final residual = 1.89681e-08, No Iterations 0
This is kOmegaSSTF
DILUPBiCG:  Solving for omega, Initial residual = 8.63384e-06, Final residual = 3.5632e-09, No Iterations 1
bounding omega, min: -39.1196 max: 2.83868e+06 average: 82692
DILUPBiCG:  Solving for k, Initial residual = 0.000213868, Final residual = 1.30345e-07, No Iterations 1
bounding k, min: -0.0417151 max: 22.2951 average: 0.703467
PIMPLE: iteration 4
DILUPBiCG:  Solving for alpha.water, Initial residual = 1.45561e-07, Final residual = 3.72721e-11, No Iterations 1
Phase-1 volume fraction = 0.991645  Min(alpha.water) = 0.0154314  Max(alpha.water) = 1
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Liquid phase volume fraction = 0.991625  Min(alpha.water) = 0.0154314  Max(alpha.water) = 1
DILUPBiCG:  Solving for Ux, Initial residual = 1.53741e-05, Final residual = 1.0046e-08, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 2.80869e-05, Final residual = 2.3609e-08, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 7.86704e-08, Final residual = 2.95063e-10, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 2.12339e-08, Final residual = 2.34074e-10, No Iterations 1
GAMGPCG:  Solving for p_rgh, Initial residual = 1.88666e-08, Final residual = 1.88666e-08, No Iterations 0
This is kOmegaSSTF
DILUPBiCG:  Solving for omega, Initial residual = 4.02313e-06, Final residual = 1.46929e-09, No Iterations 1
bounding omega, min: -22.9509 max: 2.83865e+06 average: 82692.7
DILUPBiCG:  Solving for k, Initial residual = 0.000139817, Final residual = 4.32277e-08, No Iterations 1
bounding k, min: -0.0272965 max: 22.301 average: 0.703634
PIMPLE: iteration 5
DILUPBiCG:  Solving for alpha.water, Initial residual = 1.44821e-07, Final residual = 3.73634e-11, No Iterations 1
Phase-1 volume fraction = 0.991646  Min(alpha.water) = 0.0154317  Max(alpha.water) = 1
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Liquid phase volume fraction = 0.991625  Min(alpha.water) = 0.0154317  Max(alpha.water) = 1
DILUPBiCG:  Solving for Ux, Initial residual = 3.97181e-06, Final residual = 3.51636e-09, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 9.34745e-06, Final residual = 1.30711e-08, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 5.11634e-08, Final residual = 2.30832e-10, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 2.08343e-08, Final residual = 2.0801e-10, No Iterations 1
GAMGPCG:  Solving for p_rgh, Initial residual = 1.8864e-08, Final residual = 1.8864e-08, No Iterations 0
This is kOmegaSSTF
DILUPBiCG:  Solving for omega, Initial residual = 2.51697e-06, Final residual = 1.89748e-09, No Iterations 1
bounding omega, min: -50.392 max: 2.83865e+06 average: 82692.4
DILUPBiCG:  Solving for k, Initial residual = 8.56674e-05, Final residual = 6.59347e-07, No Iterations 1
bounding k, min: -0.0407519 max: 22.3026 average: 0.703485
PIMPLE: iteration 6
DILUPBiCG:  Solving for alpha.water, Initial residual = 1.44401e-07, Final residual = 3.73939e-11, No Iterations 1
Phase-1 volume fraction = 0.991646  Min(alpha.water) = 0.0154318  Max(alpha.water) = 1
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Liquid phase volume fraction = 0.991625  Min(alpha.water) = 0.0154318  Max(alpha.water) = 1
DILUPBiCG:  Solving for Ux, Initial residual = 1.59784e-06, Final residual = 2.50779e-09, No Iterations 1
DILUPBiCG:  Solving for Uy, Initial residual = 5.68396e-06, Final residual = 3.10249e-08, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 6.23817e-08, Final residual = 2.25845e-10, No Iterations 1
GAMG:  Solving for p_rgh, Initial residual = 2.04297e-07, Final residual = 1.95536e-10, No Iterations 1
GAMGPCG:  Solving for p_rgh, Initial residual = 1.87125e-08, Final residual = 1.87125e-08, No Iterations 0
This is kOmegaSSTF
DILUPBiCG:  Solving for omega, Initial residual = 2.74728e-06, Final residual = 7.57324e-10, No Iterations 1
bounding omega, min: -28.6162 max: 2.83864e+06 average: 82692
DILUPBiCG:  Solving for k, Initial residual = 0.000144201, Final residual = 4.12669e-08, No Iterations 1
bounding k, min: -0.0272926 max: 22.3035 average: 0.703574
PIMPLE: converged in 6 iterations
ExecutionTime = 0.6 s  ClockTime = 1 s

I looked on the forum and outerCorrectorResidualControl was the new command for residualControl but OpenFOAM v-1912 does not recognize that function at all. I also checked to see if bounding k and omega are responsible for this but even their values don't change significantly. I am not able to understand why is the calculation running till 6 iterations when it clearly converges in 5 iterations itself?
For reference, here is my fvSchemes

/*--------------------------------*- 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;
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default              Euler;
}

interpolationSchemes
{
    default              linear;
}

divSchemes
{
    default              none;
    div(rhoPhi,U)        Gauss linearUpwind grad(U);
    div(phi,omega)       Gauss linearUpwind grad(omega);
    div(phi,k)           Gauss linearUpwind grad(k);
    div(phi,omega)		 Gauss linearUpwind grad(omega);
    div(phi,alpha)       Gauss vanLeer;
    div(phirb,alpha)     Gauss linear;
    div(phi,nuTilda)     Gauss linearUpwind grad(nuTilda);
    // div((muEff*dev(T(grad(U))))) Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; 
}

gradSchemes
{
    default              Gauss linear;
}

laplacianSchemes
{
    default              Gauss linear limited corrected 0.5;
}

snGradSchemes
{
    default              limited corrected 0.5;
}

fluxRequired
{
    default              none;
    p_rgh;
    pcorr;
    alpha.water;
}
wallDist
{
    method meshWave;
}

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

Please help me out here

Best
TurbulentGuy

[doerschlott]

Hello TurbulentGuy,

this might arise from your under-relaxation. You put the factor 0.8 for basically all quantities, which is then used between the PIMPLE correctors. However, as it is not specifically stated, the last corrector will not use under-relaxation. So to speak, "UFinal", "kFinal", ... are set to 1.

The solver probably even recognises the time step is converged after 5 iterations and calculates one last corrector step (#6) with no under-relaxation.


This is just a guess and worth a try, if it still appears with no under-relaxation at all. I have to admit, that I am also confused about the "residualControl" and "outerResidualControl" stuff.

One last thing for the turbulence, as I have similar problems right now. The bounding of omega looks bad to be honest. If the values differ that much, there is probably something wrong going on in the calculated physics. Just as a side note.

Good luck and let me know, if that was the source of error.

Best regards
Doerschlott

[maike151]

Hi TurbulentGuy,


I am not sure if it's the same with PIMPLE, but for SimpleFoam and in OpenFoam version 7 the syntax is like this



SIMPLE
{
nNonOrthogonalCorrectors 10;
nCorrectors 5,
nOuterCorrectors 5;
consistent yes;
residualControl
{
p 1e-4;
U 1e-4;
}
}




You can see if it worked right after you started the solver. It'll say



"Create time

Create mesh for time = 0


SIMPLE: Convergence criteria found
p: tolerance 0.0001
U: tolerance 0.0001


......"



Best
Maike

[Carlo_P]

Hey, on which version are you running?
In the openfoam7, it was a bug.


In opefoam 8 is working.


Next trick...
the residual control should be written like this:


residualControl
{
p 1e-4;
U 1e-4;
".*" 1e-4;

}


so you are sure that you will include all the variables

[TurbulentGuy]

Originally Posted by doerschlott

Hello TurbulentGuy,

this might arise from your under-relaxation. You put the factor 0.8 for basically all quantities, which is then used between the PIMPLE correctors. However, as it is not specifically stated, the last corrector will not use under-relaxation. So to speak, "UFinal", "kFinal", ... are set to 1.

The solver probably even recognises the time step is converged after 5 iterations and calculates one last corrector step (#6) with no under-relaxation.


This is just a guess and worth a try, if it still appears with no under-relaxation at all. I have to admit, that I am also confused about the "residualControl" and "outerResidualControl" stuff.

One last thing for the turbulence, as I have similar problems right now. The bounding of omega looks bad to be honest. If the values differ that much, there is probably something wrong going on in the calculated physics. Just as a side note.

Good luck and let me know, if that was the source of error.

Best regards
Doerschlott

Hello Doerschlott
First of all thank you for your insights. I did not know the final underelaxation terms.
However, even after I write them, the calculation continues even after the residual control limits are reached.
I thought maybe it's because of k and omega bounding limits(as you mentioned) but they don't change much per iteration. Do you have an idea how to take care of that?

Best
TurbulentGuy

[TurbulentGuy]

Originally Posted by Carlo_P

Hey, on which version are you running?
In the openfoam7, it was a bug.


In opefoam 8 is working.


Next trick...
the residual control should be written like this:


residualControl
{
p 1e-4;
U 1e-4;
".*" 1e-4;

}


so you are sure that you will include all the variables

Hi maike151 and Carlo_P
I am using OpenFOAM-v1912. And it does work as you guys say(the instructions are printed when I run the program).

Best
TurbulentGuy

[doerschlott]

Hey TurbulentGuy,


the thing is, that as soon as the convergence criteria is met, a final corrector loop (iteration #6) is conducted, regardless of the stated under-relaxation. Though I'm glad I could teach you about the final terms.
The residual control definitely worked, as you put 100 for the outerCorrectors and it clearly finishes the time step before that.


For the bounding problem:
I had the problem when I increased the time steps too much. I was calculating a steady-state with a transient solver and as soon as the flow field is developed, you can go well beyond CFL>1. But if you overstretch, things start to oscillate like the turbulence fields, in which case OpenFOAM starts to bound the solutions for them. As far as I know, this happens only to turbulence fields and not for "normal" fields like velocity, temperature and so on.
As the CFL number is in the range of ~1 in your case, I can only assume the boundary conditions of turbulence fields as the problem. You might check and compare them with tutorial cases of SST. Although I never simulated multi-phase with OpenFOAM, so my knowledge is very limited there. Sorry.


Best regards
doerschlott

[TurbulentGuy]

Hello doerschlott
Thank you for helping to clear my confusion.
Regarding the bounding problem, I agree with you. After your comment, I ran a small calculation on adjustable deltaT which the calculation showed around 5e-6. I am currently running the calculation at deltaT of 1e-6 and I do not see any bounding limits. I also think that having a high bound limit is fine to avoid a divide-by-zero situation.
Please correct me if I am wrong.

Warm regards
TurbulentGuy

[doerschlott]

Hey TurbulentGuy, sorry for the late response
This is true to some degree, although the real problem is having omega at near zero-values, as you divide by omega. Have a look into "kOmegaSSTBase.C", function "correctNut":

// Correct the turbulence viscosity

this->nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));

Best regards
doerschlott