CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (https://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Calculation of integral boundary values (https://www.cfd-online.com/Forums/openfoam-solving/58109-calculation-integral-boundary-values.html)

Joern Beilke (Beilke) January 12, 2005 04:01

How can I calculate the -
 
How can I calculate the

- mass flux averaged pressure for a boundary region
- mass flow rate for a region (sum, inflow, outflow)

during the calculation?

Eugene de Villiers (Eugene) January 12, 2005 06:54

Have a look at this code examp
 
Have a look at this code example:

/********************************************/
scalar volumeFluxWeightedPressureInlet = 0.0;
scalar volumeFluxInlet = 0.0;
word inletPatchName = "inlet";

//access boundary elements
const fvPatchList& patches = mesh.boundary();

//loop over boundaries
forAll(patches, patchI)
{
<BLOCKQUOTE> const fvPatch cPatch = patches[patchI];

//check boundary name
if(cPatch.name() == inletPatchName)
{
<BLOCKQUOTE> //Access boundary face area vectors
vectorField inletAreaVectors = cPatch.faceAreas();

//loop over all faces of selected
//boundary
forAll(cPatch, faceI)
{
<BLOCKQUOTE> //calculate boundary face flux
scalar cFaceFlux
= U.boundaryField()[patchI][faceI] & inletAreaVectors[faceI];

//sum face fluxes
volumeFluxInlet += cFaceFlux;

//sum flux weighted face pressures
volumeFluxWeightedPressureInlet += p.boundaryField()[patchI][faceI] * cFaceFlux;</BLOCKQUOTE></BLOCKQUOTE><BLOCKQUOTE> }</BLOCKQUOTE><BLOCKQUOTE> volumeFluxWeightedPressureInlet \= volumeFluxInlet;</BLOCKQUOTE></BLOCKQUOTE><BLOCKQUOTE> }</BLOCKQUOTE>}
/**********************************************/


Here "mesh" is the fvMesh grid, U is the velocity field (vector), p is the pressure field and "&" is the dot product operator.

The same general procedure can be used to access, modify and otherwise manipulate all types of boundary field values. For example, it can also be used to specify inlet velocities through face centre locations and or time dependent functions.

See "src/OpenFOAM/lnInclude/fvPatch.H" and its base classes in the OpenFOAM directory or the Doxygen documentation for more relevant access functions.

Eugene

Joern Beilke (Beilke) January 12, 2005 15:32

Hi Eugene, I tried to get
 
Hi Eugene,

I tried to get it running but I still get one error message:

error: 'const class Foam::fvPatch' has no member named 'faceAreas'

I think that I got an idea what you are doing in the example, but I'm a total c++ beginner ;-) I will also send you the example via email.

Thanks very much for your help

Jörn

Eugene de Villiers (Eugene) January 12, 2005 19:22

Try the function Sf() instead
 
Try the function Sf() instead (faceAreas is a polyPatch member function, which I mistakenly thought would be inherited by fvPatch).

For a C++ beginner, the Doxygen documentation can be extremely valuable; providing an easy to understand summary of class member functions and inheritance. It also has a built-in search engine, which saved me loads of time when I first started to use Foam and was still struggling with the basics of C++.

In response to your e-mail, yes "rho.boundaryField()[patchI][faceI]" will give you the density at a specified face on the specified boundary patch.

Eugene

Henry Weller (Henry) January 13, 2005 04:50

I don't agree with Eugene's e
 
I don't agree with Eugene's example as a method of calculating the patch mass flow rate. It is the flux field phi which contains the mass flux for compressible codes or the volumetic flux for incompressible codes. Also it is not usually neccessary to loop over faces; there is a large set of field functions included in FOAM which sould be used in preference to low-level looping where ever possible. Thus the inlet mass flux can simply be obtained from:

scalar massFlux = sum(phi.boundaryField()[inletPatchi]);

where the inletPatchi is the index of the inlet patch which can be obtined by searching as in Eugene's example or from

label inletPatchi = mesh.boundaryMesh().findPatchID(nameOfInlet);

where nameOfInlet is the name of the inlet patch which could be provided from a dictionary etc.

Joern Beilke (Beilke) January 13, 2005 08:33

Thanks Eugene and Henry, a
 
Thanks Eugene and Henry,

after removing some typos and changing to the Sf() function the compilation and calculation work perfectly well.

I already had a look into the Doxygen documentation but for a c++ beginner it looks a bit like "chinese backwards" ;-)

For those who want to use it I will include the running version here:

/********************************************/

scalar massFluxWeightedPressureInlet = 0.0 ;
scalar massFluxInlet = 0.0 ;
word inletPatchName = "INLET";

//access boundary elements
const fvPatchList& patches = mesh.boundary();

//loop over boundaries
forAll(patches, patchI)
{

const fvPatch& cPatch = patches[patchI];

//check boundary name
if(cPatch.name() == inletPatchName)
{

//Access boundary face area vectors
const vectorField& inletAreaVectors = cPatch.Sf();

//loop over all faces of selected
//boundary
forAll(cPatch, faceI)
{

//calculate boundary face mass flux
scalar cFaceFlux
= U.boundaryField()[patchI][faceI] & inletAreaVectors[faceI] *
rho.boundaryField()[patchI][faceI];

//sum face mass fluxes
massFluxInlet += cFaceFlux;

//sum mass flux weighted face pressures
massFluxWeightedPressureInlet += p.boundaryField()[patchI][faceI] * cFaceFlux;

}

massFluxWeightedPressureInlet /= massFluxInlet;

}

}

Info<< "MassFlowRate: " << massFluxInlet << " kg/s" << endl;
Info<< "Pressure : " << massFluxWeightedPressureInlet << " Pa " << endl;

/**********************************************/

Henry Weller (Henry) January 13, 2005 08:46

I still think this is an abso
 
I still think this is an absolutely terrible way of doing things in FOAM. I posted a much more elegant and appropriate way of calculating the patch-flux and I must warn you again that for outlets in particular the flux obtained from the velocity in this way does not obey continuity.

Joern Beilke (Beilke) January 13, 2005 09:28

Ok Henry, I replaced the
 
Ok Henry,

I replaced the

//calculate boundary face mass flux
scalar cFaceFlux
= U.boundaryField()[patchI][faceI] & inletAreaVectors[faceI] *
rho.boundaryField()[patchI][faceI];

by

scalar cFaceFlux
= phi.boundaryField()[patchI][faceI];

Your suggestion works for calculating the flux without looping the faces. But is there a way to calculate the mass flow averaged pressure without this loop?

PS: the pressure difference in now below 2% with ke and ud.

Henry Weller (Henry) January 13, 2005 10:15

You can do nearly anything in
 
You can do nearly anything in FOAM without looping. Imagine you have a field psi and want a flux-weighted average of it over patch patchi:

sum(phi.boundaryField()[patchi]*psi.boundaryField()[patchi])/sum(phi.boundaryField()[patchi])

Nils Smeds (Smeds) January 13, 2005 11:08

I guess that "...nearly anyth
 
I guess that "...nearly anything in FOAM without looping" is supposed to be read as "...nearly anything in FOAM without explicit looping"

Jarrod Sinclair (Sinclair) March 3, 2005 23:46

I am using the following code
 
I am using the following code to summarise the flux at each boundary patch, along with the net, after each timestep.

{
// loop through all boundary patches
int patchItr;
label patchIndex;
scalar flux, netFlux = 0.0;
forAll (mesh.boundaryMesh(), patchItr)
{

// get the patch index
patchIndex = mesh.boundaryMesh()[patchItr].index();

// compute the mass flux and net flux
flux = sum(phi.boundaryField()[patchIndex]);
netFlux += flux;

// compute and display the mass flux with patch name
Info<< "Mass flux at "
<< mesh.boundaryMesh()[patchIndex].name()
<< " = " << flux
<< endl;
}
Info<< "Net mass flux = " << netFlux << endl << endl;
}

This works fine, but in parallel it can't resolve some of the boundary names that are not available on the master processor (giving a name of "procBoundary0to1" for example). How can the master processor get access to these names via the mesh.boundaryMesh()[patchIndex].name() statement? (I couldn't find a solution in the Doxygen notes). However, the flux appears to be calculated fine.

Henry Weller (Henry) March 4, 2005 03:39

"procBoundary0to1" is the pro
 
"procBoundary0to1" is the processor-processor "boundary" patch which isn't a true boundary, it's the internal patch created by the decomposition and probably not important to your flux calculation. In order to sum over the "real" boundary patches in parallel you will need to test that the patches you are interested are on the processor and after the loop do a reduce operation to sum up the contributions from the processors to the boundary regions. Have a look through src/cfdTools/adjustPhi/adjustPhi.C which uses the logic you need to sum inlet and outlet fluxes to correct mass-balance for steady-state cases.

Jarrod Sinclair (Sinclair) March 7, 2005 03:11

Thanks for that Henry, I have
 
Thanks for that Henry, I have made those modifications and am using the reduce command to sum the contributions from each processor. Is there a way to get a list of all boundaries on the original mesh while running in parallel? Or is this not possible due to running on decomposed grids in seperate processes? At the moment I am defining a patch list explicity, but I could potentially parse the original constant/polyMesh/boundary file to automate this (assuming the files are available during runtime).

Henry Weller (Henry) March 7, 2005 04:48

You can send the information
 
You can send the information about the boundary patches to the master, sort the lists and send them back. However, the patch numbering on the processors will be different because each only hold the list of it's own patches so I am not sure how this complete list will help you.

sinclair March 8, 2005 22:02

I think I've got a solution th
 
I think I've got a solution that works in both serial and parallel, and uses the communcation and sorting you suggested.

Firstly, you need to #include "SortableList.H" as one of the top level includes in a solver application.

Then #include "buildGlobalBoundaryList.H" once "createMesh.H" and other initialisations have been run in the main function. (This only needs to be run once.) The file is attached below.

http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif buildGlobalBoundaryList.H

The result of this is a wordList variable called "globalPatchNames" that is available on each processor, and contains a sorted list of all real patches in the model.

To compute mass fluxes after each timestep, #include "computeMassFlux.H" (see attachment).

http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif computeMassFlux.H

I've tested this in both serial and parallel, so hopefully it works in other people's cases.

henry March 9, 2005 04:51

Thanks for these pieces of sam
 
Thanks for these pieces of sample code, it all seems fine except for one point, did you need to include Pstream.H? If so it should be included outside the "main" function at the top of the code although it is usually already included via some other includes.

sinclair March 16, 2005 20:07

You are right about Pstream.H
 
You are right about Pstream.H Henry, I've made some changes to the above code. buildGlobalBoundaryList.H is now a function that returns a wordList of patch names rather than inline code, so to avoid any variable redefinition problems. This also makes the headers a bit cleaner. However, when I perform a sort on a SortableList with OpenFOAM 1.1, I get a segmentation fault (this worked in 1.0.2). Do you know why this is so?

henry March 17, 2005 04:31

The difference between 1.0.2 a
 
The difference between 1.0.2 and 1.1 with respect to sortable list is that I have changed all calls to the C-library qsort and replaced them the equivalent STL sort calls which is certainly a more elegant method but we have had problems with it in the past.

Could you please send us the code that produces the failure and we will look into it.

sinclair March 17, 2005 09:08

Sure, the new version of "buil
 
Sure, the new version of "buildGlobalBoundaryList.H" is attached below:

http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif buildGlobalBoundaryList.H

This should be included outside the main function after top level includes like "fvCFD.H". It can then be called inside the main function like so,

wordList globalBoundaryList;
globalBoundaryList = buildGlobalBoundaryList(mesh);

Then, mass fluxes can be calculated via the following include (similar to the previous example).

http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif computeMassFlux.H

I've commented out the sort() operation in "buildGlobalBoundaryList.H" for now.

mattijs March 17, 2005 09:58

Hello Jarrod, we located fo
 
Hello Jarrod,

we located found the problem in SortableList. It is because setSize when it grows the list does not adapt the list of indices. Until we fix it you can work around it by creating your list big enough and shrinking it (i.e. setSize with something smaller than the current size)

Mattijs

mattijs March 17, 2005 10:10

This is the fix. http://ww
 
This is the fix.

http://www.cfd-online.com/OpenFOAM_D...s/mime_txt.gif SortableList.C

ali April 19, 2005 13:25

The algorithm posted by Jarrod
 
The algorithm posted by Jarrod a few threads above (for evaluating flux at the patches such as inlet etc) works fine. But what if I have a channel of say length L=1 (from x=0 to x=1) and I want to calculate flux through the middle point (x=0.5). How do I introduce a face at x=0.5 which is not a boundary patch? I found the velocities at the middle to be more accurate, so I thought calculating flux at the middle would be better. inlet flux changes its sign while it shouldn't be like that. I guess something may be wrong there.

cedric_duprat January 14, 2009 12:29

Dear all, Keeping the idea
 
Dear all,

Keeping the idea of Ali, I need to calculate a mass flow average throught a plane in my geometry.
I'll like also to define this plane (coordinate, ect...) after my calculation because I can't guess the result.
I use these few lines to create a plane inside the geometry where I can get every value I need.
// Plane 1
point pnt(0,0,0);
vector direction(0,0,1);
plane pl1(pnt,direction);
cuttingPlane cutPlane1(mesh,pl1);
const labelList& cutCells1 = cutPlane1.cells();
word setName("someCells");
cellSet currentSet1(mesh, setName, cutCells1);

// Create mesh subsetting engine
fvMeshSubset subsetter1(mesh);
label patchI = -1;
subsetter1.setLargeCellSubset(currentSet1, patchI, true);

But, how can I use them to define a patch like boundary patch which wouldn't be at the boundary of my geometry ?
So I could use previous piece of code to get different value like mass flow average or pressure average, ect ...

Thank you for your help,

Cedric

mabinty September 9, 2009 08:44

dear all!!

added the posted codes to chtMultiRegionFoam to get the net mass flux in the fluid region(s). it compiles but giving me the following warnings:

************************************************** *******************************
fluid/setRegionFluidFields.H: In function ‘int main(int, char**)’:
fluid/setRegionFluidFields.H:1: warning: unused variable ‘mesh’
fluid/setRegionFluidFields.H:5: warning: unused variable ‘K’
fluid/setRegionFluidFields.H:6: warning: unused variable ‘U’
fluid/setRegionFluidFields.H:7: warning: unused variable ‘phi’
fluid/setRegionFluidFields.H:8: warning: unused variable ‘g’
fluid/setRegionFluidFields.H:10: warning: unused variable ‘turb’
fluid/setRegionFluidFields.H:11: warning: unused variable ‘DpDt’
fluid/setRegionFluidFields.H:14: warning: unused variable ‘psi’
fluid/setRegionFluidFields.H:15: warning: unused variable ‘h
************************************************** *******************************

but the code is running. should I stil care about the warnings??

greatly appreciate your comments!

aram

sankarv April 21, 2010 15:10

Is this mass flux or mass flow rate ?
 
Hello,

I have a question about the code you have posted here.
Does it calculate Mass Flux or Mass flow rate ?

Thanks
Vaidya


Quote:

Originally Posted by mattijs (Post 202068)
Hello Jarrod,



we located found the problem in SortableList. It is because setSize when it grows the list does not adapt the list of indices. Until we fix it you can work around it by creating your list big enough and shrinking it (i.e. setSize with something smaller than the current size)

Mattijs


immortality March 30, 2013 15:11

hi
I have a question about about where should we add the includes and their related headers?

haakon March 31, 2013 07:12

immortality: Have you heard about the phrase thread necromancy before?

I really think it is bad practice to post the same question twice on the same forum at the same time. Post it once, and wait for answer.

immortality March 31, 2013 08:06

i didn't received a related answer yet.

immortality March 31, 2013 08:12

i didn't received a related answer yet.


All times are GMT -4. The time now is 11:45.