CFD Online Discussion Forums - interFoam with centrifugal Force

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)

- - interFoam with centrifugal Force (https://www.cfd-online.com/Forums/openfoam-solving/104116-interfoam-centrifugal-force.html)

interFoam with centrifugal Force

Dear Foamers,i try to simulate the development of a liquid parabola in a cylinder which is filled with water by one third. The other part is captured with air. For regarding the rotation of the fluid i included the centrifugal force in the solver interFoam like it is shown below.

Quote:

fvVectorMatrix UEqn (
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
+ rho*(omega_ ^ (omega_ ^ (mesh.C() + rad_) )) //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U)) ) & mesh.Sf())) );

For the first simulation i wanted to verify the rotational effects on the fluid. Therefore i set rad_ to zero and the angular velocity to (0 0 13.206). The gemoetry was a simple cylinder with a height of 0,5 m and a diameter of 0,15m. The liquid height is 0,2 m. After running the case and getting a constant liquid profile, i compared the solution with the analytical one and the solution calculated with fluent. While Fluent caputred the analytical solution quite well, i got an deviation of more than 2 mm in OF. So my question is where this discrepancies come from?Thx a lot!

Hello MatP,

I think the formula in your Ueqn() is wrong. Why did you add rad_ ? How did you define rad_ and omega_?

Actually I just have a little difference between MRFInterFoam and interFoam including centrifugal force for steady state and moderate angular velocities.

Modified UEqn:

Code:

    

fvVectorMatrix UEqn

    (

        fvm::ddt(rho, U)

      + fvm::div(rhoPhi, U)

      - fvm::laplacian(muEff, U)

      - (fvc::grad(U) & fvc::grad(muEff))

      //+ Fcent

      + (rho*((2.0*M_PI/60.0*omega_vector)^((2.0*M_PI/60.0*omega_vector)^mesh.C()))) //omega_vector[1/s] in [rad/m]

    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))

    );

To define omega_vector add this part to createFields.H:

Code:

dimensionedVector omega_vector(twoPhaseProperties.lookup("omega_vector"));

After compiling the solver, the magnitude and direction of the omega_vector can be set in constant/transportProperties like:

Code:

omega_vector omega_vector [ 0 0 -1 0 0 0 0 ] (0 0 5);

I hope this helps

Thanks a lot for your info styleworker, it was so helpful for my case. I had to add the centrifugal force to MRFInterFoam.

Quote:

Originally Posted by styleworker (Post 403823)

Actually I just have a little difference between MRFInterFoam and interFoam including centrifugal force for steady state and moderate angular velocities.

Modified UEqn:

Code:

    

fvVectorMatrix UEqn

    (

        fvm::ddt(rho, U)

      + fvm::div(rhoPhi, U)

      - fvm::laplacian(muEff, U)

      - (fvc::grad(U) & fvc::grad(muEff))

      //+ Fcent

      + (rho*((2.0*M_PI/60.0*omega_vector)^((2.0*M_PI/60.0*omega_vector)^mesh.C()))) //omega_vector[1/s] in [rad/m]

    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))

    );

To define omega_vector add this part to createFields.H:

Code:

dimensionedVector omega_vector(twoPhaseProperties.lookup("omega_vector"));

After compiling the solver, the magnitude and direction of the omega_vector can be set in constant/transportProperties like:

Code:

omega_vector omega_vector [ 0 0 -1 0 0 0 0 ] (0 0 5);

I hope this helps

Hi,styleworker:
I tried your method, but the code can't be successfully compiled.

Code:

omega^(omega^mesh.C())

the omega was defined

Code:

omega omega [ 0 0 -1 0 0 0 0 ] ( 0 0 0.01335 );

Code:

dimensionedVector omega

    (

        transportProperties.lookup("omega")

    );

What's missing in this code? the compile returns error 1

Without any error message, we can't give you an advice.

If you have just inserted this peace of code in your equation, your dimensions probably don't match. The dimension of the in the UEqn should be [kg/m^2s^2], because it is a volume force [N/m^3].

Your volume force should be implemented like this:

Code:

    fvVectorMatrix UEqn

    (

        fvm::ddt(rho, U)

      + fvm::div(rhoPhi, U)

      - fvm::laplacian(muEff, U)

      - (fvc::grad(U) & fvc::grad(muEff))

      + (rho*(omega ^ (omega ^ mesh.C()))) //Centrifugal Acceleration

    );

Quote:

Originally Posted by styleworker (Post 479944)

Code:

    fvVectorMatrix UEqn

    (

        fvm::ddt(rho, U)

      + fvm::div(rhoPhi, U)

      - fvm::laplacian(muEff, U)

      - (fvc::grad(U) & fvc::grad(muEff))

      + (rho*(omega ^ (omega ^ mesh.C()))) //Centrifugal Acceleration

    );

Hi,styleworker:
Thank you, but my UEqn is without rho, just looks like fvm::ddt(U), so this is no the exact reason. I have solved this problem by identify the centrifugal force as a vector field, so this works well.

I've got a private message concerning this topic. In more recent releases of interFoam transportProperties isn't initiated explicitly in createFields.H.

Therefore this code snippet has to be implemented before dimensionedVector omega is set:

Code:

IOdictionary transportProperties

(

     IOobject

     (

          "transportProperties",

          runTime.constant(),

          mesh,

          IOobject::MUST_READ,

          IOobject::NO_WRITE

     )

);

Another good tutorial is: https://openfoamwiki.net/index.php/H...Make_a_Tornado