CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   strange behavior in interPhaseChangeFoam (https://www.cfd-online.com/Forums/openfoam-solving/154468-strange-behavior-interphasechangefoam.html)

huangxianbei June 15, 2015 10:23
strange behavior in interPhaseChangeFoam
 
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.
http://www.cfd-online.com/Forums/mem...cavitation.jpg
The result is in acceptable agreement with exp

huangxianbei June 15, 2015 10:29
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 conditions:(the 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:
Code:
INLET
    {
        type            fixedValue;
        value          $internalField;
    }

    OUTLET
    {
        type            inletOutlet;
        inletValue      $internalField;
    }

  UPWALL
    {
        type            zeroGradient;
    }

    BOTTOMWALL
    {
        type            zeroGradient;
    }
   
    PART_1
    {
        type            zeroGradient;
    }
   
    frontAndBackPlanes
    {
        type            empty;
       
    }
epsilon:
Code:
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:
Code:
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:
Code:
    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
Code:
INLET
    {
        type            zeroGradient;
    }

    OUTLET
    {
        type            fixedValue;
        value          $internalField;
    }

    UPWALL
    {
        type            buoyantPressure;
    }

    BOTTOMWALL
    {
        type            buoyantPressure;
    }

    PART_1
    {
        type            buoyantPressure;
    }

    frontAndBackPlanes
    {
        type            empty;
       
    }
U
Code:
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 June 15, 2015 10:37
the transportProperties:
Code:
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:
Code:
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
Code:
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 May 18, 2017 04:46
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 May 18, 2017 20:46
Quote:
Originally Posted by l.eonardo (Post 649369)
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

Code:
max(p - pSat(), p0_)/max(p - pSat(), 0.01*pSat()),
to
Code:
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 May 21, 2017 15:28
Differences of InterPhaseChangeFoam between OF 2.3 and OF 4.1
 
2 Attachment(s)
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


All times are GMT -4. The time now is 00:23.