CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Minimal hacking to interFoam to add funny gravity

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By Phicau

Reply
 
LinkBack Thread Tools Display Modes
Old   March 9, 2016, 17:08
Default Minimal hacking to interFoam to add funny gravity
  #1
New Member
 
Daniel Duque
Join Date: Jan 2011
Location: ETSIN, Madrid
Posts: 16
Rep Power: 7
dduque is on a distinguished road
Hi,

I had this idea to model a Rayleigh-Taylor instability for the simplest possible case. I needed a funny gravity field that would pull up on one of the phases, and down on the other. I hacked the interFoam solver minimally, and I think I have got correct results. My first trial was actually succesful. UEqn.H was changed to:

Code:
    MRF.correctBoundaryVelocity(U);

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U)
            // fake gravity field
            + (2*alpha1-1) * rho * g
    );

    UEqn.relax();


    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
             (
                    mixture.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
//                  - fvc::snGrad(p_rgh)
                  - fvc::snGrad(p)
              ) * mesh.magSf()
            )
         );

        fvOptions.correct(U);
    }
Not much of a change, I introduce the funny gravity and make sure it's p, not p_rgh, the one whose gradient is evaluated.

But then I thought adding the gravity on the right-hand side of the equation would be more correct:

Code:
    MRF.correctBoundaryVelocity(U);

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U)
    );

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
             (
                    mixture.surfaceTensionForce()
                  - ghf*fvc::snGrad(rho)
//                  - fvc::snGrad(p_rgh)
                  - fvc::snGrad(p)
              ) * mesh.magSf()
            )
            // fake gravity field
            + (2*alpha1-1) * rho * g
         );

        fvOptions.correct(U);
    }
Well, no, with the latter version the system does not evolve at all. I wonder why. More importantly, I wonder if I am doing the right thing. Results are ok, but they look somewhat squearish when compared with a Lagrangian (SPH-like) simulation
comparison_step_60.jpg

Thanks in advance for any help!
dduque is offline   Reply With Quote

Old   March 9, 2016, 21:54
Default
  #2
Senior Member
 
Pablo Higuera
Join Date: Jan 2011
Location: Singapore
Posts: 339
Rep Power: 8
Phicau is on a distinguished road
Hi Daniel,

probably you have not switched on the momentum predictor in fvSolution, it is usually switched off by default.

Best,

Pablo
Phicau is offline   Reply With Quote

Old   March 10, 2016, 06:09
Default
  #3
New Member
 
Daniel Duque
Join Date: Jan 2011
Location: ETSIN, Madrid
Posts: 16
Rep Power: 7
dduque is on a distinguished road
Quote:
Originally Posted by Phicau View Post
Hi Daniel,

probably you have not switched on the momentum predictor in fvSolution, it is usually switched off by default.

Best,

Pablo
Thanks, Pablo.

You are right, it would seem that " if (pimple.momentumPredictor()) " would require the momentum predictor step. However, I have tried enabling it and now the velocity seems to evolve ... but not the pressure, which does not change, or the color field.

Aside from " momentumPredictor yes;" at the PIMPLE section, I had to add the following to fvSolutions.

Code:
  UFinal
    {
        $U;
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-06;
        relTol          0;
    }
What I don't get is, that equation contains the standard rho ( grad p + g ) term, so it should be solved no matter what. Right ??

[It's actually grad( p - rho g z ) + rho z grad(rho), but it's the same ]
dduque is offline   Reply With Quote

Old   March 10, 2016, 08:12
Default
  #4
Senior Member
 
Pablo Higuera
Join Date: Jan 2011
Location: Singapore
Posts: 339
Rep Power: 8
Phicau is on a distinguished road
Hi again,

not really, you should always include it in the first one (to assemble the system of equations properly), and if you want a momentum predictor, also in the second equation.

Momentum predictor just solves an intermediate step and updates U. Then it applies that value to obtain the independent terms (explicit) and solve the PPE. Otherwise, the U used is the one from the beginning of the time step.

Best,

Pablo
dduque likes this.
Phicau is offline   Reply With Quote

Old   March 10, 2016, 09:03
Default
  #5
New Member
 
Daniel Duque
Join Date: Jan 2011
Location: ETSIN, Madrid
Posts: 16
Rep Power: 7
dduque is on a distinguished road
Quote:
Originally Posted by Phicau View Post
you should always include it in the first one (to assemble the system of equations properly), and if you want a momentum predictor, also in the second equation.
Thanks, that seems to work. I find little difference with the previous results, but that may be ok. To check incompressibility, I am running foamCalc div phi, and the resulting field is between about -3e-4 and 3e-4. I guess that's small. Perhaps all this could go in another thread.
dduque is offline   Reply With Quote

Reply

Tags
interfoam;ueqn, rayleigh-taylor

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
InterFoam: Add an equation to only one phase voingiappone OpenFOAM Programming & Development 36 July 30, 2016 20:39
Add source term in alphaEqn.H of interFoam tayo OpenFOAM 1 October 23, 2013 03:40
InterFoam: add a source term in alpha eq. Alucard OpenFOAM 9 January 13, 2012 11:16
how to add source-term in momentum equation for interFoam? anon_g OpenFOAM 9 October 18, 2011 12:47
add temperature equation to interFoam nygbook OpenFOAM Running, Solving & CFD 3 September 4, 2011 07:28


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