How to calculate CoG in interDyMFoam
 

Hello,
I would like to calculate the centre of gravity of liquid fraction in a sloshing simulation.
Are there any options to calculate it in OF or in Paraview?

Regards
Angel

Hi angel.

i have never done it by myself, but i would use the threshold filter to selct just the cells with alpha_water > 0.5, and than use the calculator filter, in conjunction with each cell's coordinates and volumes, to obtain a new array with just one value.


good luck

You mean by calculating the weighted average of position of all the cells? Like this:

sum(pos()*alpha1*vol())/(vol()*alpha1)

In swak4Foam there is a functionobject swakExpression where the line

expression "sum(pos()*alpha1*vol())/(vol()*alpha1)";

will print exactly that at every timestep (OK. you've got to add a bit more to the controlDict, but it is less than 20 lines)

Hello,
First at all, many thanks to Bernhard and Vieri for your quick response.
I’ve tried first, with swak4foam and after installed and read yours presentations, I finally discovered one easy way to obtain the cog that works for me, but not very smart, because in can’t remove the floating point exception.
This is my controlDict, with the different options that I have tried, and the results obtained, included the error:
libs (
"libsimpleSwakFunctionObjects.so"
"libswakFunctionObjects.so"
);
functions
{
liquid_COG //simpleFunctionObjects:)
{
type swakExpression;
valueType internalField;
setName liquidCog;
accumulations (sum);
expression "(pos()*alpha1*vol())/0.034167"; //works
/*# Time sum
0.01 (0.000309988 0.0903394 0.25)*/
verbose true;
outputControl outputTime;
}
liquid_COG1 :)
{
type swakExpression;
valueType internalField;
setName liquidCog;
accumulations (sum);
expression "pos()"; //works
/*# Time sum
0.01 (-401.694 14571.8 22412.5)*/
verbose true;
outputControl outputTime;
}
liquid_COG2 :)
{
type swakExpression;
valueType internalField;
setName liquidCog;
accumulations (sum);
expression "alpha1*vol()"; //works
/*# Time sum
0.01 0.034167*/
verbose true;
outputControl outputTime;
}
liquid_COG3 :)
{
type swakExpression;
valueType internalField;
setName liquidCog;
accumulations (sum);
expression "pos()*alpha1*vol()"; //works
/*# Time sum
0.01 (1.05913e-05 0.00308663 0.00854175)*/
verbose true;
outputControl outputTime;
}
liquidCOG4:confused:
{
type swakExpression;
//functionObjectLibs ("libswakFunctionObjects.so"); Not needed if it is called in libs ("libsimpleSwakFunctionObjects.so" "libswakFunctionObjects.so");
valueType internalField;
setName liquid_Cog;
accumulations (sum);
//expression "(pos()*alpha1*vol())/(alpha1*vol())"; //Floating point exceptio
variables (
"posX=(pos().x>0.0001) ? pos().x : 0;"
"posT=vector(posX,pos().y,pos().z);"
);
expression "posT*alpha1*vol()/(alpha1*vol())"; //ERROR
/*Expression liquidCOG4 : #0 Foam::error::printStack(Foam::Ostream&) in "/home/anmartin/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/home/anmartin/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/lib/libOpenFOAM.so"
#2 __restore_rt at sigaction.c:0
#3 Foam::tmp<Foam::GeometricField<Foam::Vector<double >, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<Foam::Vector<d ouble>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/home/anmartin/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/lib/libfiniteVolume.so"
#4 parserField::FieldValueExpressionParser::parse() in "/home/anmartin3/OpenFOAM/anmartin3-2.0.x/platforms/linux64Gcc45DPOpt/lib/libswak4FoamParsers.so"
#5 Foam::FieldValueExpressionDriver::parse(std::strin g const&) in "/home/anmartin3/OpenFOAM/anmartin3-2.0.x/platforms/linux64Gcc45DPOpt/lib/libswak4FoamParsers.so"
#6 Foam::swakExpressionFunctionObject::write() in "/home/anmartin3/OpenFOAM/anmartin3-2.0.x/platforms/linux64Gcc45DPOpt/lib/libsimpleSwakFunctionObjects.so"
#7 Foam::simpleFunctionObject::execute(bool) in "/home/anmartin3/OpenFOAM/anmartin3-2.0.x/platforms/linux64Gcc45DPOpt/lib/libsimpleFunctionObjects.so"
#8 Foam::functionObjectList::execute(bool) in "/home/anmartin/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/lib/libOpenFOAM.so"
#9 Foam::Time::run() const in "/home/anmartin/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/lib/libOpenFOAM.so"
#10 main in "/home/anmartin/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/bin/interDyMFoam"
#11 __libc_start_main in "/lib64/libc.so.6"
#12 Foam::UOPstream::write(char) in "/home/anmartin/OpenFOAM/OpenFOAM-2.0.x/platforms/linux64Gcc45DPOpt/bin/interDyMFoam"
Floating point exception*/
verbose true;
outputControl outputTime;
}
};

By the other hand, in paraview, it is possible to obtain de volume with meshQuality, and then pass celldata to point data and obtain the cellcenters. Then is possible using the calculator to multiply coord(x y z)*alpha*vol, but I don’t know how to make the summation of all cell values.:confused:

If anybody knows how it is possible to solve these questiongs, I will be very gratefully.
Best regards,

I think you'll want to write sum(alpha1*vol()) instead of alpha1*vol(). Because if in ANY cell the alpha1 is 0 then you get a division by zero there. Apart from that sum(pos()*alpha1*vol())/sum(alpha1*vol) is NOT the same as sum(pos()*alpha1*vol()/(alpha1*vol)) (which is what you wrote)

Hi,

You can use function Filter -> Integrate Variables to make the summation.

Cheers,

Duong

Hello,

many thanks to everybody that helped me.

Work done perfectly.

Thanks again,
:)

A little utility
 

Hi,

I created a utility to do just this, basically as a first step to get to know the coding style, but it works for me; it reads the alpha fraction for each time step and provides the centre of mass as a weighted average (with the cell volume), and the total volume and outputs these numbers to stdout.

I'll leave this just here in case somebody could use it. It's as simple as the ptot utility on which it is based, and it seems to work too on dynamic meshes. No guarantees ;)

(unpack in $WM_PROJECT_USER_DIR/applications/utilities/postProcessing/miscellaneous , go to alpha1Tot directory and run wmake - tested on v2.1.0)