Perfect fluid implementation - interFoam
 

Hi all.

This is my first post here, so I'm going to introduce myself: my name's Gaetano and I'm doing a PhD in Chemical Engineering in Naples (Italy).

I've been lurking those forums for at least 6 months and I found answers to almost all my questions. This one, however, seems to be a bit more tricky.

I'm trying to have interFoam working with one of the two phases being a perfect fluid, i.e. having zero density and viscosity (I'd like to simply put "0" in the dictionary transportProperties). In doing so, I need to avoid all the division by density, if any. I found some in class twoPhasesMixture and I created a new class twoPhaseMixturePerfect with a little (tricky) workaround in it. When I tried to solve the damBreack tutorial with my interPerfectFoam I found this error at the very first step:

------------------------------------------------------------------------
Starting time loop

Courant Number mean: 0 max: 0
Interface Courant Number mean: 0 max: 0
deltaT = 0.00119048
Time = 0.00119048

MULES: Solving for alpha1
Phase-1 volume fraction = 0.130194 Min(alpha1) = 0 Max(alpha1) = 1
MULES: Solving for alpha1
Phase-1 volume fraction = 0.130194 Min(alpha1) = 0 Max(alpha1) = 1
#0 Foam::error::printStack(Foam::Ostream&) in "/share/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/share/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2 __restore_rt at sigaction.c:0
#3 Foam::divide(Foam::Field<double>&, double const&, Foam::UList<double> const&) in "/share/OpenFOAM/OpenFOAM-2.1.1/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#4 main in "/local/home/iitcrib3/GaetanoDM/damBreak/interPerfectFoam"
#5 __libc_start_main in "/lib64/libc.so.6"
#6 Foam::regIOobject::writeObject(Foam::IOstream::str eamFormat, Foam::IOstream::versionNumber, Foam::IOstream::compressionType) const in "/local/home/iitcrib3/GaetanoDM/damBreak/interPerfectFoam"
Floating point exception
------------------------------------------------------------------------

It seems that there is a problem with a division (see #3 above) in the p_rgh solver (as the next step should be "DICPCG: Solving for p_rgh ..."), but my insight into the problem ends here.

So, here's my question: what does the above error mean? Is there any division with a density in the denominator anywhere in interFoam?

But, besides the answers, I'm also looking for advice from expert foamers like many here: is there any chance that my attempt to "hack" interFoam could fit in its "architecture"? I mean: does the algorithm require a division by a density? I'd like to modify only the "shell" and not the "core" of interFoam (sorry for the poor analogy: I can't figure out any better way to make my point).

Thanks in advance,
Gaetano

Are you sure that the density should be 0?! As far as I know (and can find) the definition of a perfect fluid is only that the viscosity is 0, not that the density is 0.

I think that would be an Euler fluid, but I might have chosen a wrong name for my fluid. I confirm that I'm looking for something with no inertia and no viscosity. :D

Ok, in that case I would suggest using interTrackFoam from the OpenFOAM-extend project. If you define a single mesh for fluid 1 and no mesh for fluid 2, you are basically assuming fluid 2 to have 0 density and 0 viscosity.

And a small detail with respect to the density of the fluid. You will end up dividing by zero, since you create a term, which out of memory, is called rAU in the pEqn.H file. This is essentially an interpolation of your diagonal coefficients in your matrix for the momentum equation. These coefficients scale by rho (density).

The _inverse_ of this field is acting as a diffusion coefficient for the pressure on the Poisson equation, hence you are indeed dividing by zero.

Kind regards,

Niels

I am not too expereinced in OpenFOAM, but can you do the ol'e 0~ some tiny tiny number? That has worked for me in other modeling programs before.

-NG

While this is an idea worth trying, I suspect spurious currents to kill you here, as they scale with the density and viscosity ratio between both fluids.

Well, I've already thought about the possibility of rho=eps, but I'm wondering if it would be possible to have it *exactly* zero.

Here's where my evils turn true (bold mine):

It seems that I do have to play with the matrix coefficients, and I'd really like to avoid such a task. But still: would it be possible to rewrite the term rAU? And most important: is it worth an attempt? I'm not really a numeric guy...