CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Perfect fluid implementation - interFoam (https://www.cfd-online.com/Forums/openfoam-programming-development/112296-perfect-fluid-implementation-interfoam.html)

 Gaetano January 25, 2013 11:02

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

 michielm January 25, 2013 11:20

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.

 Gaetano January 25, 2013 11:33

Quote:
 Originally Posted by michielm (Post 404023) 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

 michielm January 25, 2013 11:44

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.

 ngj January 25, 2013 13:19

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

 NJG January 25, 2013 15:42

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

 akidess January 26, 2013 06:27

Quote:
 Originally Posted by NJG (Post 404084) 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.

 Gaetano January 26, 2013 13:54

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

Quote:
 Originally Posted by ngj (Post 404046) The _inverse_ of this field (rAU) is acting as a diffusion coefficient for the pressure on the Poisson equation, hence you are indeed dividing by zero.
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...

 All times are GMT -4. The time now is 13:09.