CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Post-Processing (http://www.cfd-online.com/Forums/openfoam-post-processing/)
-   -   finding average values for a plane in the vertical direction (http://www.cfd-online.com/Forums/openfoam-post-processing/114388-finding-average-values-plane-vertical-direction.html)

 NJG March 10, 2013 08:55

finding average values for a plane in the vertical direction

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

 gschaider March 13, 2013 19:46

Quote:
 Originally Posted by NJG (Post 412920) 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

 NJG March 15, 2013 10:16

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 March 15, 2013 11:18

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

 gschaider March 15, 2013 15:56

Quote:
 Originally Posted by NJG (Post 414236) 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)

 NJG March 15, 2013 16:02

Indeed. It is really just me shying away from Python for lack of experience. Thanks again for the input.

-NG

 tfuwa March 18, 2013 23:06

Quote:
 Originally Posted by gschaider (Post 414295) 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

 gschaider March 19, 2013 06:31

Quote:
 Originally Posted by tfuwa (Post 414844) 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)

 tfuwa March 19, 2013 10:43

Quote:
 Originally Posted by gschaider (Post 414934) 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)
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

 gschaider March 19, 2013 11:16

Quote:
 Originally Posted by tfuwa (Post 415004) 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 (Post 415004) 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)

 tfuwa March 25, 2013 03:44

Quote:
 Originally Posted by gschaider (Post 415019) 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

 gschaider March 25, 2013 06:55

Quote:
 Originally Posted by tfuwa (Post 416128) 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 (Post 416128) 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)

 tfuwa March 25, 2013 11:24

Quote:
 Originally Posted by gschaider (Post 416185) 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

 gschaider March 25, 2013 12:22

Quote:
 Originally Posted by tfuwa (Post 416263) Hi Bernhard, Thanks. Just let you know the wiki-page has been updated to include this usage.
Just saw it. Great

 All times are GMT -4. The time now is 23:39.