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!