CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Modelling falling solid sphere using interFoam VOF model (http://www.cfd-online.com/Forums/openfoam-solving/99158-modelling-falling-solid-sphere-using-interfoam-vof-model.html)

eelcovv March 27, 2012 10:52

Modelling falling solid sphere using interFoam VOF model
 
1 Attachment(s)
Dear Foamers,

I have been trying to model a falling perspex sphere in water using interFoam. I figured that choosing a large surface tension and dynamic viscosity would allow to represent the perspex sphere by the second liquid phase: the high surface tension should ensure a spherical shape and the high dynamic viscosity should prevent internal fluid circulations in the sphere.

With my first attemp I find that the terminal falling velocity is way too small compared to what I expect. I assume a 3 mm sphere (density 2200 kg/m3) in water (1000 kg/m3), hence I expect to find a velocity of about 0.33 m/s However, I find a velocity of about 1 mm/s, which is an order 100 too small. Please find attached a view on the velocity field and pressure distribution over the particle (bubble).

My quesions are:
1) Is it in priciple allowed to use VOF for falling object such as small particles?
2) Has anybody simulated solid spheres with VOF ?
3) Could anybody comment on my numerical settings?

My impression is that the difficulty is the large pressure inside of the bubble due to the 2*sigma/R bubble pressure. I have choosen sigma as small as possible (0.1 N/m) in order to keep the internal bubble pressure as low as possible. Nevertheless, for R=3 mm this would still lead to a P=2*0.1/3e-3=67 N/m. The hydrostatic pressure over the bubble (which takes care of the buyancy force) is 30 Pa. The pressure I find in the bubble is actually higher than anticipated: about 400 Pa. Perhaps the large pressure drop over the bubble interface give problems ?

Well, anyway, in case anybody can say anything sensible about it. Please let me know. I will include a summary of my numerical and physical settings below.

Some remarks: I have already varried some things. I checked a different gradScheme interpolation (see fvSchemes): cellMDLimited insteat of Gauss lnear. I have varied the resolution already. This mesh already uses 0.5 mlj cells. The grid was made with blockMesh and some grid refinement in the bubble areay with snappyHexMesh. Also I have tried a large surface tenstion (1 N/m), but none of it leads to a higher falling velocity

Well. That's it. Any suggestions appreciated!

Regards
Eelco


constant/transportProperties
Code:


phase1
{
    transportModel  Newtonian;
    nu              nu [0 2 -1 0 0 0 0] 1.0e-06;
    rho            rho [1 -3 0 0 0 0 0] 1000;
  sigmaC          sigmaC  [-1 -3  3 0 0  2 0 ] 22;
}

phase2
{
    transportModel  Newtonian;
    nu              nu [0 2 -1 0 0 0 0]  1;
    rho            rho [1 -3 0 0 0 0 0] 2200;
    sigmaC          sigmaC  [-1 -3  3 0 0  2 0 ] 1e-10;
}

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

0/U
Code:

    top
    {
        type            slip;
    }
    bottom
    {
        type            pressureInletOutletVelocity;
        value          uniform ( 0 0 0 );
    }
    front
    {
        type            slip;
    }
    back
    {
        type            slip;
    }
    ambient
    {
        type            slip;
    }
    wall
    {
        type            slip;
    }

0/p_rgh

Code:

boundaryField
{
    top
    {
        type            fixedValue;
        value          uniform 0;
    }
    bottom
    {
        type            buoyantPressure;
    }
    wall
    {
        type            buoyantPressure;
    }
    back
    {
        type            buoyantPressure;
    }
    ambient
    {
        type            buoyantPressure;
    }
    front
    {
        type            buoyantPressure;
    }
}

0/alpha1
Code:

boundaryField
{
    top
    {
        type            zeroGradient;
    }
    bottom
    {
        type            zeroGradient;
    }
    wall
    {
        type            zeroGradient;
    }
    back
    {
        type            zeroGradient;
    }
    ambient
    {
        type            zeroGradient;
    }
    front
    {
        type            zeroGradient;
    }
}

system/controlDict
Code:

application    interFoam;

startFrom      startTime;

startTime      0;

stopAt          endTime;

endTime        100.0;

deltaT          1e-5;

writeControl    adjustableRunTime;

writeInterval  1;

purgeWrite      0;

writeFormat    ascii;

writePrecision  6;

writeCompression uncompressed;

timeFormat      general;

timePrecision  6;

runTimeModifiable yes;

system/fvSolution
Code:

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

solvers
{
    pcorr
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance      1e-10;
            relTol          0;
            smoother        DICGaussSeidel;
            nPreSweeps      0;
            nPostSweeps    2;
            nFinestSweeps  2;
            cacheAgglomeration false;
            nCellsInCoarsestLevel 10;
            agglomerator    faceAreaPair;
            mergeLevels    1;
        }
        tolerance      1e-10;
        relTol          0;
        maxIter        100;
    }

    "(p_rgh|Psi)"
    {
        solver          GAMG;
        tolerance      1e-07;
        relTol          0.01;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels    1;
    }

    "(p_rgh|Psi)Final"
    {
        $p_rgh;
        relTol          0;
    }

    "(Bf|U|T|k|epsilon|omega|R|omega)"
    {
        solver          GAMG;
        tolerance      1e-07;
        relTol          0.1;
        smoother        GaussSeidel;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 10;
        agglomerator    faceAreaPair;
        mergeLevels    1;
    }

    "(U|k)Final"
    {
        $U;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor no;
    nCorrectors    3;
    nNonOrthogonalCorrectors 2;
    nAlphaCorr      1;
    nAlphaSubCycles 4;
    cAlpha          1;
}


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

system/fvSchemes
Code:

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

ddtSchemes
{
//    default CrankNicholson 1;
    default        Euler;
}

gradSchemes
{
//    default    Gauss linear;
          default    cellMDLimited Gauss linear 1.0;
      div(rho*phi,U) Gauss linearUpwind cellLimited Gauss linear 1;

}

divSchemes
{
    div(rho*phi,U)  Gauss limitedLinearV 1;
    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss interfaceCompression;
}

laplacianSchemes
{
    default        Gauss linear corrected;
    laplacian(gamma,Psi) Gauss harmonic uncorrected;
}

interpolationSchemes
{
    default        linear;
//    gamma            Gauss harmonic uncorrected;
}

snGradSchemes
{
    default        corrected;
}

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


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


olivierG March 28, 2012 07:54

hello,

You should try with a smaller viscosity ratio, 1e6 is too much and may give youstrong parasitic current/numerical artefacts ...
so try a a viscosity ratio of 1e3. (i.e nu2 ~1e-3)

regards,
olivier

eelcovv April 2, 2012 08:10

viscosity dependence
 
2 Attachment(s)
Hi Olivier

Thanks for the remark. It indeed appears that the ratio nu_fluid/nu_bubble was choosen too large. I have run a few cases with varying bubble viscosity. Ideally the bubble viscosity is as large as possible (to mimic a solid sphere). For lower viscosities you can anticipate a different drag coefficient of the bubble, given by

Cd=(16/Re)*((1+(3mu_p)/(2*mu_f))/(1+mu_p/mu_f))

(analytical solution for Rep<1 by Hadammard,1911)

For mu_p<<mu_f this give Cd=16/Re (stokes flow of gas bubble) and for mu_p>>mu_f this gives Cd=24/Re (stokes flow for particle). Clearly, I am not in Stokes regime, nevertheless, if I use this relation I would say that for my choise of nu_p=1e-3 m2/s -> mu_p=rho*nu_p=2.2 Pa s I am in the limit of spherical particles as Cd -> 24/Re.

Well, I run a few values of nu_p; attached the graph of the position X_p and velocity U_p of each value, incl the analytical solution for a true spherical particle. As you can see, now indeed I am in the right ball park. My previous value of nu_p (of 0.01 m2/s) gave way too low terminal falling speed, but as soon you go below 2e-3 m2/s for the kinematic viscity (i.e. take nu_p/nu_f < 1000), the terminal rising speed is reasonably well predicted.

My only concern now is that the terminal falling velocity is very sensitive to the exact choise of nu_p, whereas the values taken should all be in the limit that Cd-> 24/Re, hence still lead to the same terminal falling velocity. Perhaps this is due to differences in the deformation of the sphere. I will check the influence of the surface tension. But anybody with further suggestions: any comments appreciated:-)

Regards
Eelco

sharonyue November 28, 2012 04:57

Quote:

Originally Posted by eelcovv (Post 352665)
It indeed appears that the ratio nu_fluid/nu_bubble was choosen too large. I have run a few cases with varying bubble viscosity. Ideally the bubble viscosity is as large as possible (to mimic a solid sphere).

Quote:

Originally Posted by eelcovv (Post 352665)
You should try with a smaller viscosity ratio, 1e6 is too much and may give youstrong parasitic current/numerical artefacts ...
so try a a viscosity ratio of 1e3. (i.e nu2 ~1e-3)

Hi, if you want to simulate a perspex sphere droping in water, I dont know why you can change the nu of water, afaik,the nu of water is a constant if you dont considerate T or something.

by the way, you sphere is so small, have you ever tried a larger diameter? such as 1cm steel ball?

duongquaphim November 28, 2012 06:06

Hi Eelco,

There is a group in Sweden performing quite a lot of simulations on settling of solid particles using VOF. Here is one of their paper:

A novelmultiphase DNS approach for handling solid particles in a rarefied gas

H. Ströma, b, http://cdn.els-cdn.com/sd/entities/REcor.gif, http://cdn.els-cdn.com/sd/entities/REemail.gif, S. Sasicc, http://cdn.els-cdn.com/sd/entities/REemail.gif, B. Anderssona, b, http://cdn.els-cdn.com/sd/entities/REemail.gif
http://dx.doi.org/10.1016/j.ijmultip...ow.2011.03.011

Cheers,

Duong

emirust January 16, 2013 03:37

Hey Eelcov,

You mention you perform a grid refinement in the region close to the bubble with snapphyHexMesh. How exactly are you doing this?

Thanks!

Edit: more infor can be found : http://www.cfd-online.com/Forums/ope...plet-fall.html


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