[nicasch]

Dear OpenFoamers,

after some time spent in unsuccessfully trying to implement a relatively simple thing (like gravity) into a relatively simple application (like icoFoam), I decided to post this message.
I want to simulate laminar 2D fluid column at rest in a closed cavity, in order to obtain a hydrostatic pressure distribution with no velocity. I use icoFoam with the standard cavity tutorial. I changed the application and included gravity force into momentum equation. I also changed things for the PISO to account for the additional source in calculating pressure. The code looks like this:

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

solve(UEqn == -fvc::grad(p) + gC); // gC is volVectorField (0 -9.81 0)

// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
volScalarField rUA = 1.0/UEqn.A();

U = rUA*UEqn.H();

phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi)
+ (fvc::interpolate(rUA*gC) & mesh.Sf()); // gC added to interpolated UEqn.H()

adjustPhi(phi, U, p);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}

# include "continuityErrs.H"

U -= rUA*(fvc::grad(p));
U.correctBoundaryConditions();
}

The application compiles, but after running on cavity-case I always obtain very strange results of 2 kind:

1. BC: velocity fixedValue = (0 0 0), pressure zeroGradient at all walls. After a couple of time steps I obtain a message of continuity errors ?!

--> FOAM FATAL ERROR : Continuity error cannot be removed by adjusting the outflow.
Please check the velocity boundary conditions and/or run potentialFoam to initialise the outflow.
Specified mass inflow : 3.23118e-06
Specified mass outflow : 3.23187e-06
Difference : 6.90589e-10
Adjustable mass outflow : 0

2. When I change pressure BC on the top Wall to totalPresure, I obtain a suspiciously looking pressure distribution and, more interesting, there is some flow within closed cavity of zero velocity at walls?!



I wonder if anyone would know why this happens? I believe that the implementation of gravity is OK, or not?

Best Regards.

[nicasch]

here are pressure and velocity distributions after 5 seconds

pressure



velocity

[hjasak]

No, not that easy. The closest simple solver I can show you is boussinesqBuoyantFoam:

http://openfoam-extend.svn.sourceforge.net/viewvc/openfoam-extend/trunk/Core/Ope nFOAM-1.4.1-dev/applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussine sqBuoyantFoam.C?view=markup


If you start from there and get rid of the energy equation, you will see what happens.

Basically, if you add gravity in incompressible flows in this way, you only added body force and there is no density gradient to counter-act the effect. For example, if you have a closed box with gravity, why would there be flow at all?

Hrv

[nicasch]

Hmmm... Thank You Hrv.
Well yes, of course, there should not be any flow in such situation, that is what bothers me.
I disabled TEqn in boussinesqBuoyantFoam, but it still behaves like icoFoam without gravity. P and U are zero overall.
What I need is a simple static fluid with hydrostatic pressure distribution. I saw in my email-client that I didnot receive the whole message I've posted, there are some parts missing. If you look at the post, the things in my modified icoFoam are (I believe) the same as in boussinesqBuoyantFoam, except that I added the g-term into UEgn and into interpolated UEqn.H() needed by the PISO.
It must be possible to obtain a fine hydrostatic pressure distribution in a closed cavity for incompressible fluid without density gradient.
I must miss something really bad (

[Schag]

Hello nicasch,

I wonder if you found the solution to implement such a solver. I'm at the time trying to do quite the same thing with turbFoam, but I need gravity to simulate a flow through a dam.
I don't want to use a multiphase solver because I want to minimise the control volume (and by the way the mesh...).

So, did you found something?

Best regards

[rama13]

Don't know if it can help, but I found it easier to start from interFoam and taking away all dependencies on alpha