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/)
-   -   interFoam with centrifugal Force (http://www.cfd-online.com/Forums/openfoam-solving/104116-interfoam-centrifugal-force.html)

MatP July 3, 2012 09:23

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!

styleworker January 24, 2013 06:11

Hello MatP,

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

styleworker January 24, 2013 12:30

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

edalmau March 30, 2013 14:14

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

huangxianbei March 13, 2014 22:39

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

styleworker March 14, 2014 03:54

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


huangxianbei March 14, 2014 07:48

Quote:

Originally Posted by styleworker (Post 479944)
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
    );


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.

styleworker March 25, 2015 12:16

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


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