CFD Online Discussion Forums - Cavitation around NACA hydrofoil using interPhaseChangeFoam

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - Cavitation around NACA hydrofoil using interPhaseChangeFoam (https://www.cfd-online.com/Forums/openfoam-solving/149667-cavitation-around-naca-hydrofoil-using-interphasechangefoam.html)

Thanks.
Of course, but I want to make single phase simulation of the flow with multiphase condition, to compute the aditional vorticity which is introduced by the cavitation under the same condition.

Dear All,

I am retrieving similar strange behaviour using InterPhaseChangeFoam as it is posted in this treatment. And I would like to know if some one can help me or give some tips to solve it.

The case what I am running is a 3D geometry of a Globe Valve in a pressure drive condition. I have run the case with following pressure drops:

- No Cavitation Condition (According to experimental results)
DP = 0.7, 0.8 and 1.0 bar (Pin = 1.62, 1.74 and 1.98 bar)

- Cavitation Condition (According to experimental results)
DP = 1.5 and 1.8 bar (Pin = 2.1, 2.5, 3 bar)

The cavitation model using here is Schneer-Sauer, n = 1e+14, dNuc = 1e-05, Cc =1; Cv =1

The problem is that in all of the cases there are several cells with negative pressure (fluctuating around 0), what is no physical when absolute pressure at inlet and outlet is used. In any case I get less value than saturation pressure (3500 Pa), even when I did not expect cavitation. Normally, the negative pressure (or pressure below the saturation pressure) is located at the ring intersection between pipe and the globe valve, as is shown in this figure:

https://www.dropbox.com/s/gwywsu4pjr...eDp07.png?dl=0

https://www.dropbox.com/s/gwywsu4pjr...eDp07.png?dl=0

Firstly, I was using the interPhaseChangeFoma of OF-2.3.0 which give a high negative pressure, of order of -200000 Pa Then, changing some lines in the code, I am retrieving better results but still negative ( order of -2500 Pa) as you can see in the next figure:

https://www.dropbox.com/s/jqdq8j0juk...-Pmin.png?dl=0

https://www.dropbox.com/s/jqdq8j0juk...-Pmin.png?dl=0
The changes in code are:

1. phaseChangeTwoPhaseMixture.C in line 105

Code:

Foam::tmp<Foam::volScalarField>

Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::pCoeff

(

    const volScalarField& p

) const

{

    volScalarField limitedAlpha1(min(max(alpha1_, scalar(0)), scalar(1)));

    volScalarField rho

    (

        limitedAlpha1*rho1() + (scalar(1) - limitedAlpha1)*rho2()

    );



    return

        (3*rho1()*rho2())*sqrt(2/(3*rho1()))

       *rRb(limitedAlpha1)/(rho*sqrt(mag(p - pSat()) + 0.001*mag(pSat()))); //** 0.01*pSat-> 0.001*mag(pSta) 

}

I did not expect any effect in the pressure with this change

2. alphaEqnSubCycle.H in line 31

Code:

    else

    {

        #include "alphaEqn.H"

    }



    alpha1.max(dimensionedScalar("zero", alpha1.dimensions(), 0.0));   //**

    alpha1.min(dimensionedScalar("zero", alpha1.dimensions(), 1.0));   //**



    rho == alpha1*rho1 + alpha2*rho2;

}

This change is only to bound the alpha1 between 0 and 1.

3. UEqn.H

Code:

 /* VERSION 2.1.1*/   if (pimple.momentumPredictor())

    {

        solve

        (

            UEqn

         ==

            fvc::reconstruct

            (

                fvc::interpolate(rho)*(g & mesh.Sf())

              + (

                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)

                  - fvc::snGrad(p) 

                ) * mesh.magSf()

            )

        );

// -------------------------------------------------------------------------//

/*  VERSION OF2.3.x   if (pimple.momentumPredictor())

    {

        solve

        (

            UEqn

         ==

            fvc::reconstruct

            (

                (

                    interface.surfaceTensionForce()

                  - ghf*fvc::snGrad(rho)

                  - fvc::snGrad(p_rgh)

                ) * mesh.magSf()

            )

        );

    } */

    }

This change is only to conserve the same estructure that in 2.1.1, which use p instead of p_rgh in the snGrad, what I think can affect in positive way the pressure behaviour.

The set-up of the case is:

BOUNDARY CONDITIONS

p_rgh

Inlet → TotalPressure
Outlet → fixedPressure
Wall → fixedFluxPressure

U

Inlet → pressureInletVelocity
Oultet → zeroGradient
Wall → fixedFluxPressure

Turbulence Model

kOmegaSST

ControlDict
Adaptative time step → Co = 0.9 (PISO mode with nCorr = 4)
maximum yPlus < 3

I attach the system file for more information. Attachment 41829

I will be very greatful if some can help me in order to obtain more physical results using this solver.

Kind Regards,

Hi,
Can you post some log-output?

Hi Philipp,

Thank you for replying me.

Here you have time steps of my log. There is extra information (as min max pressure) from the implementation of solver introduced above:

Code:



Courant Number mean: 0.002765 max: 0.898738

deltaT = 8.6489e-07

Time = 0.3482777



qt = 1 Restart: no CurTim= 24310 startN= 25

DILUPBiCG:  Solving for alpha.phase1, Initial residual = 5.02842e-05, Final residual = 6.1601e-08, No Iterations 1

Phase-1 volume fraction = 0.999962  Min(alpha1) = 0.013006  Max(alpha1) = 1

MULES: Correcting alpha.phase1

Liquid phase volume fraction = 0.999962  Min(alpha1) = 0.013006  Max(alpha1) = 1

DILUPBiCG:  Solving for Ux, Initial residual = 2.6906e-05, Final residual = 4.255e-09, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.000145688, Final residual = 5.84061e-09, No Iterations 2

DILUPBiCG:  Solving for Uz, Initial residual = 3.83835e-05, Final residual = 6.46781e-11, No Iterations 2

Max vDotcP: -0    Min vDotcP: -1577.04

Max vDotvP: 1682.15    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 7.81377e-06, Final residual = 4.49767e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1574.63

Max vDotvP: 1682.24    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 6.24766e-07, Final residual = 3.15494e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1574.57

Max vDotvP: 1682.43    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 2.94952e-07, Final residual = 3.28958e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1574.57

Max vDotvP: 1682.53    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 1.44687e-07, Final residual = 1.87927e-09, No Iterations 1

qt = 1 Restart: no CurTim= 24311 startN= 25

smoothSolver:  Solving for omega, Initial residual = 3.9347e-06, Final residual = 9.11554e-10, No Iterations 2

smoothSolver:  Solving for k, Initial residual = 4.26664e-05, Final residual = 6.30417e-10, No Iterations 3

bounding k, min: 5.69094e-16 max: 33.7742 average: 1.50846

Max pressure: 249424

Min pressure: -35086.3

Max velocity: 35.3398

ExecutionTime = 582869 s  ClockTime = 583548 s



 MassFlows:   OUT.1 = 0.00550569  IN = -0.00543226

fieldMinMax minmaxdomain output:

    min(p_rgh) = -35222.6 at position (-0.0345468 0.0091031 -0.0643302) on processor 2

    max(p_rgh) = 248788 at position (-0.0343493 0.00930316 -0.0649426) on processor 2

    min(p) = -35086.3 at position (-0.0345468 0.0091031 -0.0643302) on processor 2

    max(p) = 249424 at position (-0.0343575 0.00931293 -0.0649496) on processor 2



forceCoeffs forceCoeffs_object output:

    Cm    = 0.94257

    Cd    = 54.1206

    Cl    = -23.8589

    Cl(f) = -10.9869

    Cl(r) = -12.872



Courant Number mean: 0.002765 max: 0.898441

deltaT = 8.6489e-07

Time = 0.3482786



qt = 1 Restart: no CurTim= 24312 startN= 25

DILUPBiCG:  Solving for alpha.phase1, Initial residual = 5.01624e-05, Final residual = 6.05526e-08, No Iterations 1

Phase-1 volume fraction = 0.999962  Min(alpha1) = 0.0129547  Max(alpha1) = 1

MULES: Correcting alpha.phase1

Liquid phase volume fraction = 0.999962  Min(alpha1) = 0.0129547  Max(alpha1) = 1

DILUPBiCG:  Solving for Ux, Initial residual = 2.69632e-05, Final residual = 4.20864e-09, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.000145753, Final residual = 4.02108e-09, No Iterations 2

DILUPBiCG:  Solving for Uz, Initial residual = 3.84138e-05, Final residual = 6.53723e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1574.56

Max vDotvP: 1681.62    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 7.78914e-06, Final residual = 5.28696e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1572.17

Max vDotvP: 1682.33    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 6.16447e-07, Final residual = 4.98827e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1572.08

Max vDotvP: 1682.6    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 2.88262e-07, Final residual = 3.09586e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1572.08

Max vDotvP: 1682.72    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 1.36217e-07, Final residual = 1.77707e-09, No Iterations 1

qt = 1 Restart: no CurTim= 24313 startN= 25

smoothSolver:  Solving for omega, Initial residual = 3.93651e-06, Final residual = 9.08068e-10, No Iterations 2

smoothSolver:  Solving for k, Initial residual = 4.26737e-05, Final residual = 6.35238e-10, No Iterations 3

bounding k, min: 5.62527e-16 max: 33.7998 average: 1.50844

Max pressure: 269391

Min pressure: -33755.8

Max velocity: 35.345

ExecutionTime = 582909 s  ClockTime = 583589 s



 MassFlows:   OUT.1 = 0.00550573  IN = -0.00543224

fieldMinMax minmaxdomain output:

    min(p_rgh) = -33895.6 at position (-0.0343758 0.0097322 -0.0639272) on processor 2

    max(p_rgh) = 268795 at position (-0.0331139 0.0126315 -0.0609156) on processor 2

    min(p) = -33755.8 at position (-0.0343758 0.0097322 -0.0639272) on processor 2

    max(p) = 269391 at position (-0.0331139 0.0126315 -0.0609156) on processor 2



forceCoeffs forceCoeffs_object output:

    Cm    = 0.935114

    Cd    = 54.0699

    Cl    = -23.716

    Cl(f) = -10.9229

    Cl(r) = -12.7931



Courant Number mean: 0.002765 max: 0.898713

deltaT = 8.6489e-07

Time = 0.3482795



qt = 1 Restart: no CurTim= 24314 startN= 25

DILUPBiCG:  Solving for alpha.phase1, Initial residual = 5.00192e-05, Final residual = 6.0447e-08, No Iterations 1

Phase-1 volume fraction = 0.999962  Min(alpha1) = 0.0128896  Max(alpha1) = 1

MULES: Correcting alpha.phase1

Liquid phase volume fraction = 0.999962  Min(alpha1) = 0.0128896  Max(alpha1) = 1

DILUPBiCG:  Solving for Ux, Initial residual = 2.70395e-05, Final residual = 5.07551e-09, No Iterations 1

DILUPBiCG:  Solving for Uy, Initial residual = 0.000145814, Final residual = 2.27001e-10, No Iterations 3

DILUPBiCG:  Solving for Uz, Initial residual = 3.84436e-05, Final residual = 1.07789e-10, No Iterations 2

Max vDotcP: -0    Min vDotcP: -1572.04

Max vDotvP: 1681.8    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 7.74807e-06, Final residual = 8.27649e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1569.69

Max vDotvP: 1683.2    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 6.10419e-07, Final residual = 3.94837e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1569.55

Max vDotvP: 1683.5    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 2.88158e-07, Final residual = 3.24262e-09, No Iterations 1

Max vDotcP: -0    Min vDotcP: -1569.55

Max vDotvP: 1683.6    Min vDotvP: 0

GAMGPCG:  Solving for p_rgh, Initial residual = 1.38128e-07, Final residual = 2.63841e-09, No Iterations 1

qt = 1 Restart: no CurTim= 24315 startN= 25

smoothSolver:  Solving for omega, Initial residual = 3.93721e-06, Final residual = 9.06762e-10, No Iterations 2

smoothSolver:  Solving for k, Initial residual = 4.26732e-05, Final residual = 6.32552e-10, No Iterations 3

bounding k, min: 5.56502e-16 max: 33.8248 average: 1.50841

Max pressure: 245957

Min pressure: -33841.4

Max velocity: 35.3512

ExecutionTime = 582950 s  ClockTime = 583629 s



 MassFlows:   OUT.1 = 0.00550577  IN = -0.00543221

fieldMinMax minmaxdomain output:

    min(p_rgh) = -33981.1 at position (-0.0343758 0.0097322 -0.0639272) on processor 2

    max(p_rgh) = 245322 at position (-0.0341428 0.0100374 -0.0648713) on processor 2

    min(p) = -33841.4 at position (-0.0343758 0.0097322 -0.0639272) on processor 2

    max(p) = 245957 at position (-0.0341505 0.0100462 -0.0648775) on processor 2



forceCoeffs forceCoeffs_object output:

    Cm    = 0.942619

    Cd    = 54.2101

    Cl    = -23.8661

    Cl(f) = -10.9904

    Cl(r) = -12.8757

Hi Daniel. I don't see anything unusual right now.
Some ideas:
1) Please read henry's first comment about limiting of gradient schemes:
http://www.openfoam.org/mantisbt/view.php?id=1410#c3246
So you should not do that in your fvSchemes
2) Did you test the solver with some safe settings, such as first order upwind divergence schemes and so on?
3) Did you run your loop with much more pressure cycles? Now you run 4, but with let's say 8? Even your log doesn't look suspicious, I once had a problem due to too high residuals... You will probably also need to reduce tolerance of the solver then.
4) I don't know the solver you use... do you know 100% that this is actually absolute pressure?

Thank you Philipp,

I try three different set-up according with your comments:
1) Using the scheme grad(p) Gauss linear;
2) Using upwin in div schemes
3) Using 8 pressure correctors in PISO

Here you have a figure representing the minimum pressure in this three cases.

https://www.dropbox.com/s/9vc02ca43p...ssure.png?dl=0

https://www.dropbox.com/s/9vc02ca43p...ssure.png?dl=0

All of them give similar minimum pressure and always negative. I am realy confusing what can be happen.

Trying to understand what can be happen, I would like to ask you if the turbulence model can have influence on this behaviour. I am using kOmegaSST with y+ ~ 1. Using nut = nutLowReWallFunction; k =kLowReWallFunction; omega = omegaWallFunction, to use kOmegaSST in lowRe mode. Also I try with nut = nutUSpaldingWallFunction as is indicated in this bug on the official OpenCFD bug tracker http://www.openfoam.org/mantisbt/view.php?id=179#c351.

By the way, I know that you develop a turbulence model kOmegaSST in lowRe model. Did you validated? Do you think that is apropiate to use in my case?

Thank you in advance.

Daniel,

First of all: As I understand it your solver it is uncompressional. So it doesn't matter if the pressure goes below zero. Am I right?

Secondly: I just implemented the Fluent version of low-Re SST into openFoam. Yes, I think I get the same results as in Fluent. I use it for all low-Re cases.

non cavitating simulation using InterPhaseChangeFoam

Hello all,

I've been simulating a cavitation model on Naca66 airfoil using InterPhaseChangeFoam solver, and I'm pretty new to this. I was getting abrupt pressure fluctuations over my internal field and was getting very unsatisfactory Cp values for the simulation. So I wanted to check if my model is correct or not. For this i thought of using the same solver and same setup except the cavitation model. I was using the SchnerrSauer cavitation model and set all its coefficients to zero.

But it didn't work. Can anyone help me through this ?

Daniel,

I know this post has been a few years old. Just wondering what was your conclusion and solution in the end?

Thanks,
Rdf