[Maik]

Dear everyone,
the title might be a little provocative and it might be more likely that I made a mistake in the boundary settings (that's why I choose this forum and not the bug section).

My setup is a fish passage in which one segment should be made "endless" by a cyclic boundary condition. This segment has a slope, so inlet and outlet are not on the same altitude, which I assumed would automatically generate a continuous flow as soon as inlet and outlet are connected by a cyclic boundary (because given by the slope the water level is higher at the outlet).

The problem is that after about 40s the water is almost gone!

I checked if the simulation is also loosing water when setting cyclic to a wall, but that looks fine.

I would be very thankful if an expert could have a short look at it and tell me what I did wrong to manipulate physics in such a undesired way.


Here are my boundaries:
(A link with the full project is at the end of the post.)

alpha.water before setFields

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      alpha.water;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    wall
    {
        type            zeroGradient;
    }

    atmosphere
    {
        type            inletOutlet;
        inletValue      uniform 0;
        value           uniform 0;
    }
    
    cyc_half0
    {
        type            cyclic;
        neighbourPatch  cyc_half1;
    }
    
    cyc_half1
    {
        type            cyclic;
        neighbourPatch  cyc_half0;
    }
}

U

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    wall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    atmosphere
    {
        type            pressureInletOutletVelocity;
        value           uniform (0 0 0);
    }
    
    cyc_half0
    {
        type            cyclic;
        neighbourPatch  cyc_half1;
    }
    
    cyc_half1
    {
        type            cyclic;
        neighbourPatch  cyc_half0;
    }
}

p_rgh

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p_rgh;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    wall
    {
        type            fixedFluxPressure;
        value           uniform 0;
    }

    atmosphere
    {
        type            totalPressure;
        p0              uniform 0;
        U               U;
        phi             phi;
        rho             rho;
        psi             none;
        gamma           1;
        value           uniform 0;
    }
    
    cyc_half0
    {
        type            cyclic;
        neighbourPatch  cyc_half1;
    }
    
    cyc_half1
    {
        type            cyclic;
        neighbourPatch  cyc_half0;
    }  

}

Here is the full project (<1MB):
http://sun.usf.uos.de/~maik/fish_passage.zip

just execute with:
sh Allrun.sh

p.s. I know that the mesh resolution especially at inlet, outlet and walls is too low to produce realistic values, but that shouldn't create mass defects.(?)

Best regards,
Maik

[RodriguezFatz]

Hi Maik,

something with the cyclic is going wrong. Take one timestep in paraview, such as t=4. Select mesh part "internalMesh" and volume field "alpha.water" to plot.
If you have a close look at the "upper" cyclic part, then turn it around to look at the lower cyclic part you will see that alpha differs, where it should be the same.
I will try to find the reason for that.

Philipp.

[RodriguezFatz]

Ok sorry, I was wrong... these were just interpolated values from cell center. Using the wall and patches to plot it looks ok.

[danny123]

Hello Maik,

I just filled a bug report (1477). Maybe your issue is related. In my case the source vector of the matrix equation of the pressure equation (already at pcorr level) has a wrong value in a single cell. If you resolve such an equation (provided your matrix solver succeeds to do it), all the flux calculations gets corrupted. If the flux is not correct, alpha will change too. In interFoam you do not have that problem, I did test this on the dam case in the tutorial.

You may be able to see this in your case too, if you recompile PCG to plot out the source vector, you may be able to reproduce the same error.

If you know how to manipulate the matrix solver implementation or somebody who knows, I would be grateful.

Regards,

Daniel

[RodriguezFatz]

Maik,

At least it's not the cyclic b.c.
If you set the two cyclics to wall you also lose water...
But why? I can't see any difference to the dam break tutorial, except the mesh... But that is pretty coarse in both cases.

[Maik]

Hi Philipp, Hi Daniel,
Thank you for your replies!

I used the same settings as in the damBreak tutorial, and if you also saw water disappearing in the a "just wall" scenario (actually I didn't notice that up to 40s, and the defect increases with the cyclic boundary condition) I should concentrate on the mesh.
I generated it with salome and I will have a look if I will my mesh settings are known to be incompatible with openfoam - or what will happen when I change the mesh.

Thank you and I will report my results!
Best wishes,
Maik

[Maik]

Hi!
Current status is that both setups, the cyclic boundary and the "just wall" scenario loose water. The cyclic boundary scenario does it very fast.
-There seems to be no velocity pointing outward of the volume (except for cyclic boundaries).
-checkMesh says the mesh is alright
-a run with better Courant number (0.1 instead of 1) gives the same result

Honestly, I don't know what's going on here. Did anyone experienced a similar loss of water/mass?

[RodriguezFatz]

1) Create a new mesh from scratch and try it.
2) If it doesn't work try a different meshing software.

What do you think?

[danny123]

Hi Maik,

The dam break case is set up such, that the initial condition (0 folder) is not mass conservative. The error corrects itself while integrating over time. This has nothing to do with the mesh. This is why the case looses mass much more at the begining than later. The mesh is very simple and therefore cannot be the problem.

The other possible problem, that I have sent a bug report recently for, is that pcorr is correted in a single cell within the BC by a very strange method, which leads to a very starnge initilialization. This applies only if you have set a reference pressure in fvsolutions. If you leave that out, then this problem goes away.

Henry Weller has stated that AMI has an inherent mass loss problem. I will have to investigate myself where this problem stems from. The regular interFoam has no such mass loss problem if it is properly initilalized. Maybe it is if the two AMI do not cover well. But you may check this in various threads he is refering to.

Probably opening the domain is the easiest thing to do, as he suggested.

Regards,

Daniel

[Maik]

Hi!
I reconstructed my geomatry from scratch, replacing the sharp edge that divided my two water bassins by a simple wall. However, the mass loss persists. So it looks like geometry and mesh are not the reason for the loss of mass.

Daniel, can you give more details which part of the initial conditions in the damBreak scenario is not mass conservative -- and maybe even better, how to make it mass conservative?
I didn't modify fvSolution, no pressure is that here.
And I didn't use cyclicAMI, just the normal cyclic boundaries.

Best regards and that you for your help!
Maik

[danny123]

Maik,

Dam break works with interFoam, not interDymFoam. The line "p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));" in pEqn.H is not mass conservative. correctPhi.H is also not mass conservative, but this does not matter if you problem converges well. The error may disappear over time, but you should be avare that it takes a number of p_rgh iterations to cancel out the error.

interDymFopam has a similar error, just I do not know if there are other errors. I am currently testing it. For sure, you need to state "correctPhi no" in fvsolutions.

Regards,

Daniel

[Maik]

Hi Daniel,
I didn't solve my problem.

So in principle interFoam should be mass conservative, as long as the initial/boundary conditions are correct, right?

But I dind't get what you meant with non mass conservative initial conditions. I am still using the setup from the first post and I don't see what could be wrong in the initialization.
(I am using interFoam, not interDynFoam, maybe that got mixed up.)
And what about this file: pEqn.H? Is it part of the solver?

In principle the setup is very simple: damBreak scenario as starting point, enhanced to 3D and with changed geometry: one tilted and surrounding wall, cyclic inlet and outlet on different heights and an atmosphere on top.
Does anyone know how to run it without loosing mass (water)?

[Sgs]

I have a similar problem, I'm using interFoam and initialize a drop of water in the center of a cube with cyclic top and bottom (the drop should simply fall for a long time) - but after even the first timestep I have a noticeable mass loss.

I'm getting the values of alpha (cell values, not interpolated point values) from ParaView and summing over then entire domain to evaluate the mass loss. Not sure if this is the correct way to do this calculation.

[danny123]

Hi Sgs,

This has been a while and I made progress since. I will try to summarize my thinking:

In the dam break tutorial, there is an amount of heavier liquid 1 in the lower left corner of a rectangular space filled with a lighter fluid 2. At time 0, there is supposely no movement, meaning U = 0. No movement also means that phi = 0 and phiAlpha1 = 0.

Trouble is that the momentum balance on a cell (UEqn) can only be fulfilled if the pressure gradient is unidirectional in direction of gravity (and / or surface tension). This means, the surface boundary has to be horizontal (to fulfill above criterion both for gravity and surface tension). This is not the case for the right side of the alpha1 area. If you try to solve for pEqn, you will find a pressure field that is in line with div phi = 0, but not phi = 0. If phi is not 0 then U is also not 0. So, you violate your initial condition.

To overcome this problem, I ended up by putting a wall on the right hand alpha1 area boundary, calculate the pressure field and start the simulation without the artificial wall using the pre-calculated pressure field and U = 0. This seems to me like you would do in an experiment as well since you need a force to keep the alpha1 fluid in the left corner before you release the fluid.

Summarized in short: if phi = 0 then div phi 0 and U = 0. If div phi = 0, phi is not necessarily 0, consequently U also not.

If you have a drop, you may use a similar workaround. You surround your drop left and right hand with an artificial wall and calculate your pressure field. You will have a connection to the area of alpha2 at the bottom of the drop. This does not matter since the drop cannot fall down (So U = 0). At time 0, you remove the artificial walls and you will see the drop detach from the upper wall.

Regards,

Daniel