CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Track vortices using Swak4FOAM

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   April 11, 2012, 11:16
Default Track vortices using Swak4FOAM
  #1
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Hi all,

I am simulating the propagation in space and time of some vortices.
I would like to be able to trace the center of these vortices into my field and then be able to compute the circulation on a closed surface around this center.

So, first step I need to trace the center of my vortex. It looks like defining some functions in the controlDict and using the libsimpleSwakFunctionObjects could be an option but I am not sure how to do it.
To start easy I was planning to find the center of my vortex by looking for the maximum vorticity in my field.

I wrote those lines in my dictionary:
Code:
functions 
(
    defineFollowVortex {
        type addGlobalVariable;
        outputControl timeStep;
        outputInterval 1;

        globalScope followVortex;
        globalVariables {
            yThres {
                valueType scalar;
                value 0.2;
            }
        }
    };
    createVortexCenter {
        type createSampledSet;
        outputControl timeStep;
        outputInterval 1;
        setName vortexCenter;
        set {
            type cloud;
            axis y;
            points( (0.5 2 0.0) );
        }
    }

    vorticityMax {
	type swakExpression;
	verbose true;
	valueType set;
	setName vorticityMax;
	expression "(pos.y()>yThres) ? magVorticity : 0"; 
	accumulations ( max );
        set {
            type swakRegistryProxy;
            axis y;
            setName vortexCenter;
        }
        interpolate true;
        interpolationType cellPoint;
    }

    vortexHeight {
        type swakExpression;
        verbose true;
        valueType set;
        setName  fillingHeight;
        globalScopes ( "fillDam" );
        expression "(magVorticity >= 0.99*vorticityMax) ? pos().y : 0";
        accumulations (
            max
        );
        set {
            type swakRegistryProxy;
            axis y;
            setName vortexCenter;
        }
        interpolate true;
        interpolationType cellPoint;
    }
);
However I am really not sure it is the way to go and I get the following error message:
Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

[0] 
[0] 
[0] --> FOAM FATAL IO ERROR: 
[0] "ill defined primitiveEntry starting at keyword 'functions' on line 60 and ending at line 124"
[0] 
[0] file: /CFDBB/vincent/OpenFOAM/vincent-2.1.0/run/tutorials/compressible/rhoPimpleDyMFoamWithSources/LES/AfterMeeting120330/testVortexHeight/processor0/../system/controlDict at line 124.
[0] 
[0]     From function primitiveEntry::readEntry(const dictionary&, Istream&)
[0]     in file lnInclude/IOerror.C at line 132.
[0] 
FOAM parallel run exiting
Any idea?
vinz is offline   Reply With Quote

Old   April 11, 2012, 15:53
Default
  #2
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Hi all,

I am simulating the propagation in space and time of some vortices.
I would like to be able to trace the center of these vortices into my field and then be able to compute the circulation on a closed surface around this center.

So, first step I need to trace the center of my vortex. It looks like defining some functions in the controlDict and using the libsimpleSwakFunctionObjects could be an option but I am not sure how to do it.
To start easy I was planning to find the center of my vortex by looking for the maximum vorticity in my field.

I wrote those lines in my dictionary:
Code:
functions 
(
    defineFollowVortex {
        type addGlobalVariable;
        outputControl timeStep;
        outputInterval 1;

        globalScope followVortex;
        globalVariables {
            yThres {
                valueType scalar;
                value 0.2;
            }
        }
    };
    createVortexCenter {
        type createSampledSet;
        outputControl timeStep;
        outputInterval 1;
        setName vortexCenter;
        set {
            type cloud;
            axis y;
            points( (0.5 2 0.0) );
        }
    }

    vorticityMax {
	type swakExpression;
	verbose true;
	valueType set;
	setName vorticityMax;
	expression "(pos.y()>yThres) ? magVorticity : 0"; 
	accumulations ( max );
        set {
            type swakRegistryProxy;
            axis y;
            setName vortexCenter;
        }
        interpolate true;
        interpolationType cellPoint;
    }

    vortexHeight {
        type swakExpression;
        verbose true;
        valueType set;
        setName  fillingHeight;
        globalScopes ( "fillDam" );
        expression "(magVorticity >= 0.99*vorticityMax) ? pos().y : 0";
        accumulations (
            max
        );
        set {
            type swakRegistryProxy;
            axis y;
            setName vortexCenter;
        }
        interpolate true;
        interpolationType cellPoint;
    }
);
However I am really not sure it is the way to go and I get the following error message:
Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

[0] 
[0] 
[0] --> FOAM FATAL IO ERROR: 
[0] "ill defined primitiveEntry starting at keyword 'functions' on line 60 and ending at line 124"
[0] 
[0] file: /CFDBB/vincent/OpenFOAM/vincent-2.1.0/run/tutorials/compressible/rhoPimpleDyMFoamWithSources/LES/AfterMeeting120330/testVortexHeight/processor0/../system/controlDict at line 124.
[0] 
[0]     From function primitiveEntry::readEntry(const dictionary&, Istream&)
[0]     in file lnInclude/IOerror.C at line 132.
[0] 
FOAM parallel run exiting
Any idea?
I'm not sure about the syntax error. Maybe ; following defineFollowVortex.

Anyway. What might help you is minPosition/maxPosition (this is currently only available in the development-version of swak. You'll have to get that via Mercurial): they take an expression and return the position where that expression has its minimum/maximum (you'll have to guess from the names which does which)
gschaider is offline   Reply With Quote

Old   April 12, 2012, 03:08
Default
  #3
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Dear Bernhard,

Thank you for your reply. Indeed it looks like minPosition and maxPosition could be a good starting point.

I went to the wiki in order to see how to download the development version. I got it and apparently managed to install in on my openfoam version 2.1.

I looked at the examples and did not find one including the min or max Position function. Could you tell me if something like this should work:
Code:
vorticityMax {
	type swakExpression;
	verbose true;
	valueType pointField;
	setName vorticityMax;
	expression "maxPosition(magVorticity)"; 
    }
What does the function return? is it a vector with x y z coordinate?
Or do you have an example?
vinz is offline   Reply With Quote

Old   April 12, 2012, 09:08
Default
  #4
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
I got it working. At least I can track one vortex center with this code:
Code:
functions 
(
    maxVortPosY
    {
        type swakExpression;
	valueType internalField;
        expression "maxPosition(mag(vorticity)).y";
	accumulations ( 
	   max
	);
        verbose true;        
    } 
    maxVortPosX
    {
        type swakExpression;
	valueType internalField;
        expression "maxPosition(mag(vorticity)).x";
	accumulations ( 
	   max
	);
        verbose true;        
    } 

);
Now I am trying to do an XY plot instead of X versus time using pyFoamPlotWatcher.
vinz is offline   Reply With Quote

Old   April 12, 2012, 15:52
Default
  #5
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
I got it working.
Great. Sorry for a lack of an example. I tested that feature with another case (that I can't publish). But if you have a simple example I'll be happy to include it into swak

Quote:
Originally Posted by vinz View Post
At least I can track one vortex center with this code:
More will be difficult with the min/max-approach. Unless you "mask out" the region of the first vortex

Quote:
Originally Posted by vinz View Post
Now I am trying to do an XY plot instead of X versus time using pyFoamPlotWatcher.
Sorry. PyFoamPlotWatcher can currently only do plots over time. Thanks for the reminder on that missing feature (but of course you can easily plot the data written by it easily with any program you like)
gschaider is offline   Reply With Quote

Old   April 13, 2012, 01:56
Default
  #6
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Quote:
Originally Posted by gschaider View Post
More will be difficult with the min/max-approach. Unless you "mask out" the region of the first vortex
Indeed, this is my problem. I have two counter-rotating vortices propagating in my field. They move in time so they are not always at the same place
I am not sure how to handle that. I would need to say something like "track the maximum vorticity only where vorticity.z is positive/negative" but I'm not sure if there is a way to do it.

Quote:
Originally Posted by gschaider View Post
Sorry. PyFoamPlotWatcher can currently only do plots over time. Thanks for the reminder on that missing feature (but of course you can easily plot the data written by it easily with any program you like)
Nothing to be sorry of, you already do a lot! When you say that I can plot the data written by PyFoamPlotWatcher, are you talking about the file in the different directories swakExpression_*? or is there another way to use the data analysed by PyFoamPlotWatcher in order to create the graph I want?

Last edited by vinz; April 13, 2012 at 02:32.
vinz is offline   Reply With Quote

Old   April 16, 2012, 02:38
Default
  #7
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Dear Bernhard,

do you think it is possible to create a subset of my mesh based on the vorticity value or on a coordinate for instance and then look for a maximum only in that subset?
vinz is offline   Reply With Quote

Old   April 16, 2012, 04:50
Default
  #8
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Dear Bernhard,

do you think it is possible to create a subset of my mesh based on the vorticity value or on a coordinate for instance and then look for a maximum only in that subset?
That's what I meant with "mask out". Something like (this is just a sketch):

vortPos1=maxPosition(vorticity)
maskField=(mag(pos()-vortPos1)<rCut) ? 0 : 1
vortPos2=maxPosition( maskField ? vorticity : -1000)
maskField=(mag(pos()-vortPos2)<rCut) ? 0 : maskField
vortPos3= ....

Of course the criteria for masking could be something more advanced (the problem here is that rCut being to small might give you "vortex centers" that are in the other vortex and being too large might "hide" another vortex center)

Alternative would be "dynamic" cellSets based on an expression, but that feature is not yet implemented (partly because I'm not sure whether the parts of the OF-code where this could be interesting can "cope" with cellSets that change their size)
gschaider is offline   Reply With Quote

Old   April 16, 2012, 10:32
Default
  #9
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Ok I see, not straight forward but brilliant!

I manage to get something which looks interesting with the following lines:
Code:
    mask
    {
        type swakExpression;
	valueType internalField;
        variables (
            "maskField=(vorticity.z>0 || pos().y<0.3) ? 0 : 1;"
	);
        storedVariables (
            {
                name maskField;
                initialValue "0";
            }
        );
        expression "maskField";
	accumulations ( 
	   max
	);
        verbose true;        
    } 


    maxVortPos
    {
        type swakExpression;
	valueType internalField;
        variables (
	    "vortPos2= maxPosition((maskField>0) ? magVorticity : -1000);"
        );
        storedVariables (
            {
                name vortPos2;
                initialValue "vector(-0.4,1,0)";
            }
        );
        expression "vortPos2";
	accumulations ( 
	   max
	);
        verbose true;        
    } 

    maxVortPosY
    {
        type swakExpression;
	valueType internalField;
        expression "vortPos2.y";
	accumulations ( 
	   max
	);
        verbose true;        
    } 

    maxVortPosX
    {
        type swakExpression;
	valueType internalField;
        expression "vortPos2.x";
	accumulations ( 
	   max
	);
        verbose true;        
    }
Now I would like to compute the integral of U on a surface which would be a sphere centered on my vortex.
I tried the following to create such a sphere :
Code:
    distField
    {
        type expressionField;
        outputControl timeStep;
        outputInterval 1;
        fieldName distanceFromVortex;
        expression "mag(pos()-vortPos2)";
        autowrite true;        
    } 

    surfaceIsoDist
    {
        type swakExpression;
        valueType surface;
        surfaceName isoVortexDist;
        surface {
            type isoSurface;
            isoField distanceFromVortex;
            isoValue 0.15;
            interpolate true;
        }
        verbose true;
        expression "mag(U)";
        accumulations (
            min
            max
	    average
        );
    }
but it crashes complaining that vortPos2 does not exist or is of wrong type although I use it earlier to get vortPos2.x for instance.

Then I am wondering what type is expected at this place? Any idea?
vinz is offline   Reply With Quote

Old   April 16, 2012, 14:04
Default
  #10
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Ok I see, not straight forward but brilliant!
I think it is the other way round: straight forward (brute force) but not brilliant/elegant

Quote:
Originally Posted by vinz View Post

Now I would like to compute the integral of U on a surface which would be a sphere centered on my vortex.
I tried the following to create such a sphere :
Code:
    distField
    {
        type expressionField;
        outputControl timeStep;
        outputInterval 1;
        fieldName distanceFromVortex;
        expression "mag(pos()-vortPos2)";
        autowrite true;        
    } 

    surfaceIsoDist
    {
        type swakExpression;
        valueType surface;
        surfaceName isoVortexDist;
        surface {
            type isoSurface;
            isoField distanceFromVortex;
            isoValue 0.15;
            interpolate true;
        }
        verbose true;
        expression "mag(U)";
        accumulations (
            min
            max
	    average
        );
    }
but it crashes complaining that vortPos2 does not exist or is of wrong type although I use it earlier to get vortPos2.x for instance.

Then I am wondering what type is expected at this place? Any idea?
I think you misunderstood the "storedVariable": this means that the variable keeps its value between timesteps. What you want is either

a) putting vortPos2 in an expression field
b) using a global variable (see slide 35ff in http://openfoamwiki.net/images/d/da/...Leoben2011.pdf).

Something like (untested)

Code:
    
calculateVortex {
        type calculateGlobalVariables;
        outputControl timeStep;
        outputInterval 1;
        valueType internalField;
        toGlobalNamespace vortexLoc;
        toGlobalVariables (
           vortexPos2
        );
        variables (
            "vortPos2= maxPosition((maskField>0) ? magVorticity : -1000);"
        );
}
Then make the scope known with

Code:
globalScopes ( vortexLoc );
where you want to use that location. See the circulatingSplash-example
gschaider is offline   Reply With Quote

Old   April 19, 2012, 02:47
Default
  #11
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Thanks again for your help Bernhard!

I managed to get it working, and to track my two vortices.
Here is an example of what I was trying to do:
vinz is offline   Reply With Quote

Old   April 19, 2012, 07:23
Default
  #12
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Thanks again for your help Bernhard!

I managed to get it working, and to track my two vortices.
Here is an example of what I was trying to do:
Nice. Could you include this as "swak4Foam was not designed to do this - searching for vortices" on the Wiki-page?

BTW: looking at the picture it occurred to me that by searching for the minimum and the maximum of "vorticity*pos().x" the "masking" could have been avoided
gschaider is offline   Reply With Quote

Old   April 19, 2012, 10:52
Default
  #13
Senior Member
 
Vincent RIVOLA
Join Date: Mar 2009
Location: France
Posts: 277
Rep Power: 9
vinz is on a distinguished road
Honestly I have never user the wikki so I have no idea how to include this

However, it is not urgent, it's an ongoing work so I cannot share everything right now. The contract ends on july and I should have other interesting (hopefully) thing to show afterwards and place it on the wikki.

Regarding the other possibility to track the vortices, I think it is less general. And although it might work in this specific case, I have other cases with wind comming from the left or where one vortex might turn around the other and I think your last proposal would not work anymore.

Thanks again for your help.
vinz is offline   Reply With Quote

Old   April 23, 2012, 19:18
Default
  #14
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,915
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by vinz View Post
Honestly I have never user the wikki so I have no idea how to include this
That is Wiki with one K. Unless you're talking about the company

It's quite easy. Really

Quote:
Originally Posted by vinz View Post
However, it is not urgent, it's an ongoing work so I cannot share everything right now. The contract ends on july and I should have other interesting (hopefully) thing to show afterwards and place it on the wikki.

Regarding the other possibility to track the vortices, I think it is less general. And although it might work in this specific case, I have other cases with wind comming from the left or where one vortex might turn around the other and I think your last proposal would not work anymore.

Thanks again for your help.
Well. The problem with your solution is that if you're unlucky the vortices can oszillate (exchange positions): if the vorticities for both of them are almost the same then for some time vortex A has the "first maximum", then B, then A again .... and in your graph you'll see the vortices exchanging place.

But I guess there is no general solution that will always work
gschaider is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
swak4Foam for faceZones andrea.pasquali OpenFOAM 3 August 30, 2011 11:24
von Karman vortex street. Vortices formation mechanism. mnvl Main CFD Forum 1 February 24, 2010 18:53
Bluff body vortices philippose OpenFOAM Running, Solving & CFD 14 April 11, 2008 08:02
Particle track James Bewley CFX 0 April 28, 2006 06:04
Particle track files Sankalp CFX 1 August 25, 2004 23:25


All times are GMT -4. The time now is 14:05.