# Minimal hacking to interFoam to add funny gravity

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

 March 9, 2016, 17:08 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 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!

 March 9, 2016, 21:54 #2 Senior Member   Pablo Higuera Join Date: Jan 2011 Location: Singapore Posts: 339 Rep Power: 8 Hi Daniel, probably you have not switched on the momentum predictor in fvSolution, it is usually switched off by default. Best, Pablo

March 10, 2016, 06:09
#3
New Member

Daniel Duque
Join Date: Jan 2011
Posts: 16
Rep Power: 7
Quote:
 Originally Posted by Phicau 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 ]

 March 10, 2016, 08:12 #4 Senior Member   Pablo Higuera Join Date: Jan 2011 Location: Singapore Posts: 339 Rep Power: 8 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.

March 10, 2016, 09:03
#5
New Member

Daniel Duque
Join Date: Jan 2011
Posts: 16
Rep Power: 7
Quote:
 Originally Posted by Phicau 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.

 Tags interfoam;ueqn, rayleigh-taylor

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post voingiappone OpenFOAM Programming & Development 36 July 30, 2016 20:39 tayo OpenFOAM 1 October 23, 2013 03:40 Alucard OpenFOAM 9 January 13, 2012 11:16 anon_g OpenFOAM 9 October 18, 2011 12:47 nygbook OpenFOAM Running, Solving & CFD 3 September 4, 2011 07:28

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