[eelcovv]

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

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

    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

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

0/alpha1

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

system/controlDict

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

/*--------------------------------*- 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

/*--------------------------------*- 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]

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]

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]

Originally Posted by eelcovv

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).

Originally Posted by eelcovv

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]

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, , , S. Sasicc, , B. Anderssona, b,
http://dx.doi.org/10.1016/j.ijmultip...ow.2011.03.011

Cheers,

Duong

[emirust]

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

[aspsarv]

Hello ever one,

Could you please let me know how I can model a solid drop into a tank?

Which solver do you recommend?

I have three phase of solid gas and liquid. Interfoam works? what about VOF?