[Edwon]

Hi Foamers!

I am trying to do a simple water volume calculation in a two-phase flow. To find the volume of the water in inter(DyM)Foam, I can simply use

dimensionedScalar waterVolume = sum(mesh.V()*alpha1);

However, for my simulations I need to know the volume of water in a certain part of the simulation domain, for example, let's say x<0 and x>=0. One of my attempts was this:

dimensionedScalar waterVolumeLeft = sum(mesh.V()*alpha1*neg(mesh.C().internalField().component(vector::X)));
dimensionedScalar waterVolumeRight = sum(mesh.V()*alpha1*pos(mesh.C().internalField().component(vector::X)));

However, this does not result in the correct volume: as a test I add both pos(..) and neg(..), which should result in the same waterVolume, but instead it gives this in a test run:
waterVolume: 1.31441405254865e-11, left+right:7.109575809296095e-12.

I've also tried to write the x-positions to a volScalarField, and looping over the elements, but I keep getting a (consistent) volume of 7.1095...e-12, with something like this:

volScalarField xPositions
(
    IOobject
    (
        "xPositions",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    dimensionedScalar("xPositions",dimensionSet(0,1,0,0,0,0,0), 0.0)
);
forAll(xPositions,iter)
{
    xPositions[iter]=mesh.C().internalField()[iter].component(vector::X);
}

dimensionedScalar waterVolumeLeftCheck  = 0.0;
dimensionedScalar waterVolumeRightCheck = 0.0;
forAll(xPositions,iter)
{
    if (xPositions[iter]<0)
        waterVolumeLeftCheck+=mesh.V()[iter]*alpha1[iter];
    else
        waterVolumeRightCheck+=mesh.V()[iter]*alpha1[iter];
}

Can anyone point me in the right direction with this? Any help is appreciated!

Best,
Edwin

[Edwon]

Hmm after another search, I just found this topic:
http://www.cfd-online.com/Forums/ope...integrate.html

It looks like the fvc::domainIntegrate(alpha1*neg(xPositions)) function is the solution. I still don't see how my other attempts failed to give the right answer..

[chriss85]

Just a rough sketch:
scalar volume = 0;
forAll(mesh, cellI)
{
if(mesh.C().x() <= 0)
volume += mesh.V()[cellI]*alpha1[cellI];
}

[Edwon]

Originally Posted by chriss85

Just a rough sketch:
scalar volume = 0;
forAll(mesh, cellI)
{
if(mesh.C().x() <= 0)
volume += mesh.V()[cellI]*alpha1[cellI];
}

That was my first try as well, but it gives an incorrect value. I link to a solution to the problem in my previous post.

[wouter]

hello Edwon,
Did you try this with a little larger volume, because it could be a numerical accuracy problem.

hope this helps,
Wouter

[chriss85]

Did you run this on multiple cores or on a single core? The code I showed is only suitable for a single core. If you need multiple cores, you need to run reduce(volume, sumOp<scalar>()); at the end to add the results of all the separate meshes together.

[Edwon]

I just checked, that also solves the problem! (I assumed that the contributions of all cores would add up automagically)

I still find it strange that sum(mesh.V()*alpha1) works on multiple cores, and sum(mesh.V()*alpha1*neg( mesh.C().internalField().component(vector::X) )) does not. But, we (you) found out what was causing the discrepancy and that is what counts!

@Wouter, I've just checked, that was not the case. Haven't thought to look at that before, though.