CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Post-Processing

finding average values for a plane in the vertical direction

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

Reply
 
LinkBack Thread Tools Display Modes
Old   March 10, 2013, 08:55
Default finding average values for a plane in the vertical direction
  #1
NJG
Member
 
Nick Gutschow
Join Date: Jan 2013
Posts: 36
Rep Power: 4
NJG is on a distinguished road
Hi Folks,

I have a finished run that looks good, paraView is showing expected behavior, but now I wish to create some plots for my data that paraView won't do.

Specifically, I want to take some scalar, let's say T, and average it on a given horizontal plane. So let's say z is vertical and x and y go from 0 to 10; I want to average the T values for all 100 volumes (spread through my x-y plane) at a given z value. Then, I want to end up with results that I can use to plot z vs T in a conventional Cartesian manner.

When I investigate my T file for some time step however, there is just a huge list of values. A given value is not explicitly tagged with a given x,y,z value in that file.

So, my questions are as follows I guess:

1) Is anyone aware of good way to accomplish my go to produce data for my z vs T plot?

2) If not, looking ahead to making my own function, where is it referenced what index in a given scalar file at a given timestep corresponds to a given spatial location? Or, is there a trick to output or export OpenFOAM data in a more readable and manipulable format to my end?

Thanks everyone!

-NG
NJG is offline   Reply With Quote

Old   March 13, 2013, 19:46
Default
  #2
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,912
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by NJG View Post
Hi Folks,

I have a finished run that looks good, paraView is showing expected behavior, but now I wish to create some plots for my data that paraView won't do.

Specifically, I want to take some scalar, let's say T, and average it on a given horizontal plane. So let's say z is vertical and x and y go from 0 to 10; I want to average the T values for all 100 volumes (spread through my x-y plane) at a given z value. Then, I want to end up with results that I can use to plot z vs T in a conventional Cartesian manner.

When I investigate my T file for some time step however, there is just a huge list of values. A given value is not explicitly tagged with a given x,y,z value in that file.

So, my questions are as follows I guess:

1) Is anyone aware of good way to accomplish my go to produce data for my z vs T plot?

2) If not, looking ahead to making my own function, where is it referenced what index in a given scalar file at a given timestep corresponds to a given spatial location? Or, is there a trick to output or export OpenFOAM data in a more readable and manipulable format to my end?

Thanks everyone!

-NG
Well you can have the cake and eat it to: do it in paraview although paraview doesn't know it. The thing is that Paraview ... if you have a Paraview with Python support compiled in ... and the version is not totally old .... comes with numpy a nice library for matrix operations. Instead of cutting planes you'd divide the range of z into N "bins" have a look which cell falls into which bin (there is a function in numpy for that) and calculate a weighted average. Approximately like this:
Code:
bins=numpy.linspace(min(minY),max(maxY),25)
vol,edges=numpy.histogram(heightField,bins=bins,weights=volume)
alpha,edges=numpy.histogram(heightField,bins=bins,weights=volume*alphaField)

xVals=0.5*(edges[1:]+edges[0:-1])
vals=alpha/vol

data=numpy.asarray([xVals,vals]).T
numpy.savetxt(dataFile("alphaByHeight"),data)
The things why you can't immediately use this is:
- This is for alpha1 not for T
- for height Y not Z
- You'll probably have to import numpy
- This was lifted from an example of swak4Foam (fillingTheDam) and the variables minY, maxY, heightField (which is Y), volume (the cell volumes), and alphaField reach the python-namespace through some swak-Stuff (this is from a function-object that evaluates this python-code at each write-time).

I must admit I haven't done much python in Paraview recently so I can't tell you how to best get to the data (with the "Python Calculator", "Programmable Filter" or by working on the Python Shell) but if you get to a solution I think it would be much appreciated if you shared it here
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request
gschaider is offline   Reply With Quote

Old   March 15, 2013, 10:16
Default
  #3
NJG
Member
 
Nick Gutschow
Join Date: Jan 2013
Posts: 36
Rep Power: 4
NJG is on a distinguished road
gschaider,

Thanks for the reply! I have got to say however, there MUST be an easier way. Actually, this is driving me a bit nuts now.

Using the sample utility I am able to produce a 1d line of data using the "raw format for sets as follows:

1 1 1 300
2 1 1 301
3 1 1 302

Where that is x, y, z, and T for temperature say. But I can only vary one of the x, y or, z values I guess? I want an output just like this, but I want it to just pump out the T value for every cell, not just the cells in a single line.

I find it hard to believe that this functionality does not exist in this utility. It seems so basic... The most basic of outputs; just a text file giving (x y z <field 1> <field 2> <field 3>....) in columns.

Does anyone know if and how the sample utility could be used to do this?

-Nick
NJG is offline   Reply With Quote

Old   March 15, 2013, 11:18
Default
  #4
NJG
Member
 
Nick Gutschow
Join Date: Jan 2013
Posts: 36
Rep Power: 4
NJG is on a distinguished road
Ok, I think I found my own work around.

From paraView I can just save my data as a csv, which outputs all my fields and points.

From here I am half way through a MATLAB code to basically area weight average the scalar in question for a given plane and then set these to vertical displacement values equal to the center point of that plane.

-NG
NJG is offline   Reply With Quote

Old   March 15, 2013, 15:56
Default
  #5
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,912
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by NJG View Post
Ok, I think I found my own work around.

From paraView I can just save my data as a csv, which outputs all my fields and points.

From here I am half way through a MATLAB code to basically area weight average the scalar in question for a given plane and then set these to vertical displacement values equal to the center point of that plane.
Which basically boils down to my solution (or are you writing out the plane data?): binning all cells at the same height together.

That is the nice thing about people: that they have different views what is "easy" and each of these ideas is valid: for you the idea of clicking through two different programs and letting them "communicate" through files, for me to invest some work in the beginning and then just starting the run and have all the data ready for me when it is finished (of course your way is more time efficient if you only have to do this once for one run)
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request
gschaider is offline   Reply With Quote

Old   March 15, 2013, 16:02
Default
  #6
NJG
Member
 
Nick Gutschow
Join Date: Jan 2013
Posts: 36
Rep Power: 4
NJG is on a distinguished road
Indeed. It is really just me shying away from Python for lack of experience. Thanks again for the input.

-NG
NJG is offline   Reply With Quote

Old   March 18, 2013, 23:06
Default
  #7
Member
 
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 74
Blog Entries: 1
Rep Power: 6
tfuwa is on a distinguished road
Quote:
Originally Posted by gschaider View Post
Which basically boils down to my solution (or are you writing out the plane data?): binning all cells at the same height together.

That is the nice thing about people: that they have different views what is "easy" and each of these ideas is valid: for you the idea of clicking through two different programs and letting them "communicate" through files, for me to invest some work in the beginning and then just starting the run and have all the data ready for me when it is finished (of course your way is more time efficient if you only have to do this once for one run)
Hi there,

I would like to choose part of a patch (for instance, cylinder faces in a specific domain), and then calculated the forces on this selected part. Is there a way to do this in OpenFoam at run-time? Thanks a lot, Albert
tfuwa is offline   Reply With Quote

Old   March 19, 2013, 06:31
Default
  #8
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,912
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by tfuwa View Post
Hi there,

I would like to choose part of a patch (for instance, cylinder faces in a specific domain), and then calculated the forces on this selected part. Is there a way to do this in OpenFoam at run-time? Thanks a lot, Albert
I'm not clear on your definition of "cylinder faces" and "calculate forces".

If you're going for complete patches you might look into the forces-functionObject that comes with OF. If you only want to select parts of the patches then you can use for instance a swak-patch expression that looks like this:

Code:
the Forces {
  type patchExpression;
  patches (
    wall1
    wall2 
  );
  variables (
     "radius=0.1;"
     "forces= .... ;" // that is the part you have to specify
  );
  expression "(sqrt(pos().x*pos().x+pos().y*pos().y)<radius) ? forces : vector(0,0,0)";
  accumulations (
     sum
  );
}
that would sum up the force for all faces on the patches which are within 0.1 to the z-axis. Specification of force is you responsibility but you can use any field that is available and lots of operators (see the groovyBC-documentation)
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request
gschaider is offline   Reply With Quote

Old   March 19, 2013, 10:43
Default
  #9
Member
 
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 74
Blog Entries: 1
Rep Power: 6
tfuwa is on a distinguished road
Quote:
Originally Posted by gschaider View Post
I'm not clear on your definition of "cylinder faces" and "calculate forces".

If you're going for complete patches you might look into the forces-functionObject that comes with OF. If you only want to select parts of the patches then you can use for instance a swak-patch expression that looks like this:

Code:
the Forces {
  type patchExpression;
  patches (
    wall1
    wall2 
  );
  variables (
     "radius=0.1;"
     "forces= .... ;" // that is the part you have to specify
  );
  expression "(sqrt(pos().x*pos().x+pos().y*pos().y)<radius) ? forces : vector(0,0,0)";
  accumulations (
     sum
  );
}
that would sum up the force for all faces on the patches which are within 0.1 to the z-axis. Specification of force is you responsibility but you can use any field that is available and lots of operators (see the groovyBC-documentation)
Thank you, Bernhard, this is exactly what I want.

For pressure force, the expression is simple,

Code:
"forces= p*sf()*mag(sf()) ;"
But for viscous force,

Code:
"forces= rho()*nu*dev(twoSymm(grad(U))) *sf()*mag(sf());"
It seems to me that dev, grad, curl and laplacian are not available in groovyBC. How to apply those operators? Many thanks. Also, please correct me if there are any errors in the force expressions.

Regards,
Albert
tfuwa is offline   Reply With Quote

Old   March 19, 2013, 11:16
Default
  #10
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,912
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by tfuwa View Post
Thank you, Bernhard, this is exactly what I want.

For pressure force, the expression is simple,

Code:
"forces= p*sf()*mag(sf()) ;"
But for viscous force,
mag(sf()) is called area in OF. So alternative might be "p*normal()*area()" (not sure about the sign)

Quote:
Originally Posted by tfuwa View Post
Code:
"forces= rho()*nu*dev(twoSymm(grad(U))) *sf()*mag(sf());"
It seems to me that dev, grad, curl and laplacian are not available in groovyBC. How to apply those operators? Many thanks. Also, please correct me if there are any errors in the force expressions.
The reason is that to calculate these you would have to calculate them in the whole field and I don't want a lowly BC that "only" works on a patch to trigger that.

For your expression: would have to look it up myself how the original functionObjects calculate that. But what might be sufficient to reformulate it in terms of snGrad(U) (that is available in groovyBC). Otherwise you'd have to use an expressionField-functionObject that creates a field named devGradU with the expression "dev(twoSymm(grad(U)))" and use that field in your expression (not sure but you might have to use internalField(devGradU) if that field is 0 on the boundary).

Also: rho is a field. So no (). nu depends on the solver whether this is available as a field (only then swak can access it. There are ways to get it anyway)
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request
gschaider is offline   Reply With Quote

Old   March 25, 2013, 03:44
Default
  #11
Member
 
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 74
Blog Entries: 1
Rep Power: 6
tfuwa is on a distinguished road
Quote:
Originally Posted by gschaider View Post
mag(sf()) is called area in OF. So alternative might be "p*normal()*area()" (not sure about the sign)


The reason is that to calculate these you would have to calculate them in the whole field and I don't want a lowly BC that "only" works on a patch to trigger that.

For your expression: would have to look it up myself how the original functionObjects calculate that. But what might be sufficient to reformulate it in terms of snGrad(U) (that is available in groovyBC). Otherwise you'd have to use an expressionField-functionObject that creates a field named devGradU with the expression "dev(twoSymm(grad(U)))" and use that field in your expression (not sure but you might have to use internalField(devGradU) if that field is 0 on the boundary).

Also: rho is a field. So no (). nu depends on the solver whether this is available as a field (only then swak can access it. There are ways to get it anyway)
Hi Bernhard,

Many thanks for your help.

The final expression is

Code:
p*rho*normal()*area()-rho*nu*snGrad(U)*area()
where the first term is due to pressure and the second represents the contribution of viscous.

BTW, what is the difference between sf() and area() ?

Regards,
Albert
tfuwa is offline   Reply With Quote

Old   March 25, 2013, 06:55
Default
  #12
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,912
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by tfuwa View Post
Hi Bernhard,

Many thanks for your help.

The final expression is

Code:
p*rho*normal()*area()-rho*nu*snGrad(U)*area()
where the first term is due to pressure and the second represents the contribution of viscous.
Cool. Even cooler if you could add a small recipe to
http://openfoamwiki.net/index.php/Co...Usage_examples

Quote:
Originally Posted by tfuwa View Post
BTW, what is the difference between sf() and area() ?
Basically "Sf()" is "area()*normal()" (for OpenFOAM that describes all it needs to know about that face in a single value)
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request
gschaider is offline   Reply With Quote

Old   March 25, 2013, 11:24
Default
  #13
Member
 
Albert Tong
Join Date: Dec 2010
Location: Perth, WA, Australia
Posts: 74
Blog Entries: 1
Rep Power: 6
tfuwa is on a distinguished road
Quote:
Originally Posted by gschaider View Post
Cool. Even cooler if you could add a small recipe to
http://openfoamwiki.net/index.php/Co...Usage_examples



Basically "Sf()" is "area()*normal()" (for OpenFOAM that describes all it needs to know about that face in a single value)
Hi Bernhard,

Thanks. Just let you know the wiki-page has been updated to include this usage.

Kind regards,
Albert
tfuwa is offline   Reply With Quote

Old   March 25, 2013, 12:22
Default
  #14
Assistant Moderator
 
Bernhard Gschaider
Join Date: Mar 2009
Posts: 3,912
Rep Power: 40
gschaider will become famous soon enoughgschaider will become famous soon enough
Quote:
Originally Posted by tfuwa View Post
Hi Bernhard,

Thanks. Just let you know the wiki-page has been updated to include this usage.
Just saw it. Great
__________________
Note: I don't use "Friend"-feature on this forum out of principle. Ah. And by the way: I'm not on Facebook either. So don't be offended if I don't accept your invitation/friend request
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
Problem with Gmsh nishant_hull Open Source Meshers: Gmsh, Netgen, CGNS, ... 18 April 22, 2015 08:43
finding average temperature of room sunilpatil CFX 1 February 3, 2013 06:03
boundaries with gmshToFoam‏ ouafa Open Source Meshers: Gmsh, Netgen, CGNS, ... 7 May 21, 2010 12:43
Finding Average Temperature and Solid Fraction atulverma FLOW-3D 1 May 14, 2009 12:39
Finding average x-velocity at any particular x. Dpu FLUENT 2 April 11, 2006 22:32


All times are GMT -4. The time now is 08:41.