[huangxianbei]

Hi,all:
I'm doing a cavitation simulation around the hydrofoil recently with interPhaseChangeFoam. However, I found quite strange results from this solver.

Description of my model: A hydrofoil fixed in a channel , the free stream velocity (inlet velocity is U = 5.33m/s), the cavitation number is defined as : a = (p_out - p_sat)/(0.5*rho*U^2)=1.25, as in interPhaseChangeFoam, the reference pressure is set to 0 Pa by default, the pressure in the calculation is the absolute pressure. Therefore, p_out = 20055.6 Pa according to the definition of the cavitation number.

The mesh used is about 20000 elements, with yPlus ranging from 37 to 267, which is suitable to use the wall function.

First of all, I ran a case (3D, 0.5 million cells) in simpleFoam to see the pressure distribution under noncavitating condition, here is the comparison with experiment.

The result is in acceptable agreement with exp

[huangxianbei]

As 3D simulation is time consuming, I study 2D case in interPhaseChangeFoam.

Still, I'd like to see if the noncavitating case agrees with the exp. So I set the following boundary conditionsthe estimation of k and epsilon can be found http://openfoamwiki.net/index.php/TurboECPGgi2D and http://support.esi-cfd.com/esi-users/turb_parameters/)

alpha1:

INLET
    {
        type            fixedValue;
        value           $internalField;
    }

    OUTLET
    {
        type            inletOutlet;
        inletValue      $internalField;
    }

   UPWALL
    {
        type            zeroGradient;
    }

    BOTTOMWALL
    {
        type            zeroGradient;
    }
    
    PART_1
    {
        type            zeroGradient;
    }
    
     frontAndBackPlanes
    {
        type            empty;
        
    }

epsilon:

INLET
    {
        type            fixedValue;
        value           $internalField;
    }

    OUTLET
    {
        type            inletOutlet;
        inletValue      $internalField;
    }

    UPWALL
    {
        type            epsilonWallFunction;
        value           $internalField;
    }

    BOTTOMWALL
    {
        type            epsilonWallFunction;
        value           $internalField;
    }

    PART_1
    {
        type            epsilonWallFunction;
        value           $internalField;
    }

     frontAndBackPlanes
    {
        type            empty;
        
    }

k:

INLET
    {
        type            fixedValue;
        value           $internalField;
    }

    OUTLET
    {
        type            inletOutlet;
        inletValue      $internalField;
    }

    UPWALL
    {
        type            kqRWallFunction;
        value           $internalField;
    }

    BOTTOMWALL
    {
        type            kqRWallFunction;
        value           $internalField;
    }

    PART_1
    {
        type            kqRWallFunction;
        value           $internalField;
    }

     frontAndBackPlanes
    {
        type            empty;
        
    }

nut:

    INLET
    {
        type            calculated;
        value           uniform 0;
    }

    OUTLET
    {
        type            calculated;
        value           uniform 0;
    }

    UPWALL
    {
        type            nutkWallFunction;
        value           uniform 0;
    }

    BOTTOMWALL
    {
        type            nutkWallFunction;
        value           uniform 0;
    }

    PART_1
    {
        type            nutkWallFunction;
        value           uniform 0;
    }

     frontAndBackPlanes
    {
        type            empty;
        
    }

p_rgh

INLET
    {
        type            zeroGradient;
    }

    OUTLET
    {
        type            fixedValue;
        value           $internalField;
    }

    UPWALL
    {
        type            buoyantPressure;
    }

    BOTTOMWALL
    {
        type            buoyantPressure;
    }

    PART_1
    {
        type            buoyantPressure;
    }

     frontAndBackPlanes
    {
        type            empty;
        
    }

U

INLET
    {
        type            fixedValue;
        value           $internalField;
    }

    OUTLET
    {
        type            pressureInletOutletVelocity;
        phi             phi;
        value           $internalField;
    }

    

    PART_1
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    UPWALL
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    BOTTOMWALL
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
     frontAndBackPlanes
    {
        type            empty;
        
    }

[huangxianbei]

the transportProperties:

phaseChange on;
//transportModel Newtonian;
//nu nu [0 2 -1 0 0 0 0] 1e-6;
phaseChangeTwoPhaseMixture SchnerrSauer;

pSat             pSat       [1 -1 -2 0 0]    2300;   // saturation pressure

sigma           sigma [1 0 -2 0 0 0 0] 0.07;

phase1
{
    transportModel Newtonian;
    nu              nu [0 2 -1 0 0 0 0] 1e-6;
    rho             rho [1 -3 0 0 0 0 0] 1000;
}

phase2
{
    transportModel Newtonian;
    nu              nu [0 2 -1 0 0 0 0] 4.273e-04;
    rho             rho [1 -3 0 0 0 0 0] 0.02308;
}
SchnerrSauerCoeffs
{
    n               n      [0 -3 0 0 0 0 0]     1.6e+13;
    dNuc            dNuc   [0 1 0 0 0 0 0]      2.0e-06;
    Cc              Cc     [0 0 0 0 0 0 0]      1;
    Cv              Cv     [0 0 0 0 0 0 0]      1;
}

fvSchemes:

ddtSchemes
{
    default              Euler;
}

interpolationSchemes
{
    default              linear;
}

divSchemes
{
    default              none;
    div(rhoPhi,U)        Gauss linearUpwind grad(U);
    div(phi,epsilon)       Gauss linearUpwind grad(epsilon);
    div(phi,k)           Gauss linearUpwind grad(k);

    div(phi,alpha)       Gauss vanLeer;
    div(phirb,alpha)     Gauss interfaceCompression;
}

gradSchemes
{
    default              Gauss linear;
}

laplacianSchemes
{
    default              Gauss linear limited 0.5;
}

snGradSchemes
{
    default              none;
    snGrad(pd)           limited 0.5;
    snGrad(rho)          limited 0.5;
    snGrad(alpha1)       limited 0.5;
}

fluxRequired
{
    default              none;
    p_rgh;
    pcorr;
    alpha1;
}

fvSolution

solvers
{
    alpha1
    {
        maxUnboundedness 1e-5;
        CoCoeff          2;
        maxIter          5;
        nLimiterIter     2;

        solver           PBiCG;
        preconditioner   DILU;
        tolerance        1e-12;
        relTol           0.1;
    };

    "(U|k|epsilon)"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0.1;
    }

    "(U|k|epsilon)Final"
    {
        solver          PBiCG;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }

    p_rgh
    {
        solver           GAMG;
        tolerance        1e-8;
        relTol           0.1;

        smoother         DICGaussSeidel;
        nPreSweeps       0;
        nPostSweeps      2;

        cacheAgglomeration true;

        nCellsInCoarsestLevel 10;
        agglomerator     faceAreaPair;
        mergeLevels      1;

        maxIter          50;
    };

    pcorr
    {
        $p_rgh;
        relTol           0;
    };

    p_rghFinal
    {
        solver            PCG;
        preconditioner
        {
            preconditioner  GAMG;

            tolerance        1e-6;
            relTol           0;

            nVcycles         2;

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

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

PIMPLE
{
    momentumPredictor          no;
    nOuterCorrectors           1;
    nCorrectors                2;
    nNonOrthogonalCorrectors   1;

    cAlpha                     0;
    nAlphaCorr                 1;
    nAlphaSubCycles            1;
}


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

As can be seen , all the settings are the same with the tutorial: cavitatingBullet

The resulting Cp at x/l = 0.1 is about -2, which is much lower than the exp, indicating a much lower pressure than exp. While the calculation is at least stable.

I used to think the time scheme may affect the result, so I change it to the second order backward/ Crank , however, the two schemes used can easily lead to overflow>

I have got stuck on this problem for days Is there anyone who can tell me what's wrong? Any advice is appreciated.

Best regards
Xianbei

[l.eonardo]

Hallo to everyone in the Forum.

Huang,
have you figured out how to improve your model to avoid an over estimation of the cavitating phenomena? I have, more or less, the same set up you use and I obtain a cavity longer than the experimental.

Thank you in advance
Tom

[huangxianbei]

Originally Posted by l.eonardo

Hallo to everyone in the Forum.

Huang,
have you figured out how to improve your model to avoid an over estimation of the cavitating phenomena? I have, more or less, the same set up you use and I obtain a cavity longer than the experimental.

Thank you in advance
Tom

Hello:
The reference is here:
http://www.tfd.chalmers.se/~hani/kur...ChangeFoam.pdf

The major modification is:
from

max(p - pSat(), p0_)/max(p - pSat(), 0.01*pSat()),

to

max(p - pSat(),p0_)/max(p - pSat(), 0.001*mag(pSat())),

Meaning that the model is more sensitive to the vaporization or condensation.

Maybe this method can improve your results.

Xianbei

[Mashiro5]

Dear All,

I'm having similar problems in the case of simple cavitation estimation using interPhaseChangeFoam.
My test case is the well known NACA66 hydrofoil by Shen and Dimotakis (Shen YT and Dimotakis PE. The influence of surface cavitation on hydrodynamic force. In: Proceedings of the 22nd ATTC, St. Johns, NL, Canada, 8–11 August 1989, pp.44–53.).

I already simulated this hydrofoil using OpenFOAM 2.3 with quite satisfactory results.

When I moved to OpenFOAM 4.1, I re-ran some old calculations to check if everything behaves similarly.
Unfortunately this is not the case. As you can see from the figures below, using OF 4.1 leads to significantly longer cavity bubbles (longer than OF 2.3 and in turn, to the experiments).
Calculations have been carried with exactly the same setup: same mesh, same numerical schemes (except for the new conventions in fvSchemes of OF 4), same numerics (exactly the same fvSolution file). Calculations were initialized with the same non-cavitating calculations in order to provide a realistic estimation of the pressure around the hydrofoil and avoid/limit the initial transient.

Do you have any explanation circa these differences?
Many thanks,
Stefano