[paka]

As stated in topic.
I plan to use calculated pressures at each time step. The pressures which appear on the nodes of the rigid structure.

Could anyone please explain me the way and classes which are responsible for reading them out?

I would like to implement them to some numeric integration scheme, so one could calculate overall acting force on the structure.

Thanks.

[philippose]

Hello Krystian,

A Good evening to you!

I think the following utility will do if not exactly what you need...almost exactly what you need.... :-)!

It basically lets you specify a set of patches, and it calculates the the total force acting on each of those patches due to the pressure acting on those surfaces, and finally, also gives the total resultant force along the given direction vectors (along with some other data such as total surface area of each patch, and the area weighted mean pressure acting on each patch).

calcPressureForces.tar.gz

Enjoy!

Philippose

[paka]

Hi Philippose and good morning (at least here),

That's wonderful, that's exactly what I was looking for.

But... some questions.

Current calcPressureForce tool calculates forces based on the p file in time directories. I'm using two-phase flow solver (rasInterFoam) which produces pd files instead of p. I tried to dig a bit into your code, however my OpenFOAM code knowledge is limited.

Therefore, does any one know how to modify that code so it will read data from pd files? I do not know how to access it.

Now some other questions related to that tool. If anything is trivial, please forgive me, for me CFD is still a new area of research.

1.
In dictionary file (calcPressureForcesDict) there is following line:

// Density for conversion of normalised pressure to physical pressure units
// in the case of incompressible solvers
rho rho [1 -3 0 0 0 0 0] 849.19;

I assume one does not need to change anything here. This is standard conversion?
Probably, I should read some literature about it.

2.
Does it also work for 3-D cases?
If yes, I assume the only thing is to specify:

dirVector3 (0 0 1);

3.
Does this tool work similar as liftDrag tool works? I mean, to calculated lift and drag coefficients one needs to calculate forces first. Therefore, it seems that calcPressureForces is like a sub-code of liftDrag. Is it right or am I mistaken?

Thanks for all your help,
Krystian

[paka]

OK. Temporarily, I figured out some workaround. I made easy shell script, which simply makes a copy of pd files to p files.

This seems to be working, results look good. It would be nice if some OpenFOAM guru could say something if such workaround makes sense.

<u>One more question regarding to calcPressureForces</u>:
What I see there is no output file from that tool. So, if one wants to make a plot wrt time, I think I would need to modify the source code so it will produce continuous data which could be viewed by some plotting software.

Thanks,
Krystian

[paka]

Sorry for being paint in neck, but one more OpenFOAM/C++ related question.

I know people already asked couple times for writing out the output to file, but always it was kind of long-complicated for newcomers.

Maybe someone well familiar with OpenFOAM coding could provide some simple template, so anyone could use and apply one or two lines of code to any of the source files. This way, anyone could obtain additional output file the way they want.

In my case, I want to read one of the variables, in this case sum of forces on the wall, which is produced on the screen, but I just want to extract the numbers, which I can also do easily in source file. However, I want to specify just one or two lines which will do something like that:


Read $some_number - which I know how to do.

then

Write that number to $some_file.dat


That way I could obtain the file with series, like:
1
2
3
4
etc.
and then nicely plot it.

Right now I can do it using some BASH scripting, but I guess the cleaner, nicer, and faster way is to implemented it straight into any of OpenFOAM source files.

Thanks,
Krystian

[paka]

So many posts today, I hope you forgive me
But when I work, most of you sleep.

I need to know how surface area is calculated in 2D space for above tool. I do not have an access at the moment to view my model and surface force is calculated based on the area surface.

I find in the code that:
scalar totSurfArea = gSum(mesh.magSf().boundaryField()[forcePatchID]);
but cannot decipher this.

When I converted 2D model from Fluent to OpenFOAM it created somehow the "weird" thickness of around 0.0248 (or 0.248), which normally does not affect 2D computation in any way.

However, when I calculate surface forces I wonder what thickness is used? Is it a unity used for area computation, or is that 0.0248 number which is taken?

Thanks,
K

[philippose]

Hello Krystian :-)!

A Good <morning/evening/afternoon> to you (<- Choose the option right for you :-)!)

Ok... so now for the questions... shall work through them one by one...:

1. To change the name of the pressure data file from "p" to "pd", all you need to do is... go to the file "createP.H", and in the following code snip... make the given change:

IOobject
(
"p", <--- change this "p" to "pd"
runTime.timeName(),
mesh,
IOobject::MUST_READ
)


2. In the dictionary file, the factor "rho", is the density of the medium you are working with. Since I work with hydraulic oil, I have a value of "849.19", which is the density of oil. You will have to change this based on the density of the fluid you are working with.
Since you are working with an incompressible solver, just changing the value to the fluid you use, should work.... though I am not sure what happens when there are two fluids, with different densities effecting the system

3. The tool works without any modifications, for 3D cases too... you don't need to specify a "dirVector3", because the tool automatically calculates the 3rd direction to be normal to the two vectors given...

So... as the dictionary is set up currently, dirVector1 = (1 0 0), and dirVector2 = (0 1 0), which means, internally, the tool calculates the third direction to be (0 0 1).

4. Yes... this tool is in a way quite similar to the lift drag utility, but different in that, here, I limit the calculation to only the pressure forces. The fluid shear forces are not calculated.

5. The utility "calcPressureForces" puts out all the output to the screen... which is basically the output that you see when you run the utility... it does not write any files with any result data...

What I do currently to make plots, is to calculate the forces for each time interval, and save that into a file using the command:

calcPressureForces [root] [case] > datafile

Then I use the linux "grep" command to extract the data lines I need, and the linux "colrm" command to remove all the extra data, so that only the numbers of interest remain... this I finally plot with "gnuplot"

If you need the utility to directly write the data to a file, you will need to modify the code.

6. The total surface area of a patch is calculated using the command:

scalar totSurfArea = gSum(mesh.magSf().boundaryField()[forcePatchID]);

In this case....
* gSum is the command to sum up values... something similar to a "sigma" in mathematics

* mesh.magSf() gives you the surface area of each face of each cell present in the mesh.

* mesh.magSf().boundaryField()[forcePatchID] gives you the surface area of each face present in the patch specified by the patch ID "forcePatchID".

A sum of the surface area of all the faces on a given patch gives you the surface area of that patch.

7. If your patch is located such that the side with thickness 0.0248 is selected, then the thickness will be used in the calculation of the surface area, because you cannot calculate the "surface area" of a one-dimensional "face" :-)!

Hope this helps you go further with your work...!

Have a nice day!

Philippose

[paka]

Hi Philoppose!

Many very thanks for you detailed answers. They clarify a lot!

Let me provide some additional questions and comments below. I will use numbers, as before, so they will help relate comments to particular questions.

2a.
First of all, I come from the structural filed, I'm a structural engineer, so all fluid stuff is, as I probably said before, very new for me.

So, from structural point of view, having pressure calculated, to obtain forces acting on the wall structure I would simply integrate the pressure over the area on which that pressure acts and this way I would obtain forces.

How is your computation done? Why do you use rho definition? Isn't it already included in pressure? I tried to track the code, but no success.

2b.
Yes, I work with incompressible solver, but with two-phase flow solver. Right now, the only thing which I can do and be more or less sure about result is the time step where all structure is fully covered with water. Other way it is a bit insane to pick particular steps and adjust them to where water (gamma field) is.

Here comes my question.
Do you know anyway or function/classes which are responsible for gamma field? I would like to use gamma filed to adjust domain density. So calcPressureForce tool could check the gamma filed and based on gamma value assign corresponding density values.
Do you think this makes sense? That way the tool might be extended to two-phase flows.

5.
I did exactly the same way as you did. Only exception, I modified the code so it generates directly the input files for my plotting software;
btw. if someone there is a Mac user I recommend the "Plot" software, very nice, high quality 2D plots.

6,7.
Thanks for detailed explanation. I very appreciate it. Helped to understand the code. Even if I look at wiki page, somehow I get lost, maybe that is because everything is still new and also I never used C++ before.

Question.
What pressure integration method you use? As I mentioned above, I would simply integrate the pressure to obtain acting forces. In this case simple trapezoidal or rectangular rule should work without trouble.
I did some simple test on one of my small patches, however, results which I obtain from trapezoidal rule are much different.
For example, the output from calcPressureForce tool is in Newtons.. which at time step which I check gives me 95818[N] force. However if I use trapezoidal rule, the force in my case is 95[N]. So I guess I missed something.
I even thought that this is maybe the trouble between air/water densities (1/1000), however that doesn't make a sense because on the plot whole structure is fully submerged, so that cannot be a trouble.

OK. I should be finishing.. here is 11PM And I'm still in my office.

Very many thanks for all the input and help!

Best regards,
Krystian

[paka]

If I could also ask, the code is:

patchPressureForce = gSum(p.boundaryField()[forcePatchID]
* mesh.Sf().boundaryField()[forcePatchID]
);

meanPressure = gSum(p.boundaryField()[forcePatchID]
* mesh.magSf().boundaryField()[forcePatchID]
);


First of all, what is the difference between those two guys? It seems they are exactly the same. Even in output they give exactly the same values.

Second of all, maybe I'm wrong, but it seems it is:
p * A
That looks like force.. hmm... OK, time to leave from work.

[philippose]

Hello Krystian :-)!

A Good day to you! Hope you were able to sleep well, inspite of the headache of trying to understand the code :-)!

Ok... here are the answers to your queries:

2a. When you use an incompressible fluid solver, one of the assumptions made, is that the density of the fluid remains a constant at all points in time and space within the domain of interest.

So.... for the incompressible solvers in OpenFOAM that I use, the pressure specified, is not the actual pressure (for example in Pa), but normalised pressure. The pressure is normalised with density. The unit of "normalised" pressure is:

Pa/rho = (N/m²)/(kg/m³) = [kg/(ms²)]/(kg/m³) = (m²/s²)

Hence, if you need to specify a pressure of say.... 1 bar at some boundary, you would have to convert it using:

1 bar => (1 * 1.013e+5)/(density)

And that converted value would be the value you would have to use for the boundary condition in the solver (in my case, simpleFoam, icoFoam, etc..single phase solvers).

This means.... the pressure field in the file "p", is not the real pressure, but the normalised pressure. Now, if you need to calculate the force due to the calculated pressure on a surface, you need to multiply the calculated "normalised" pressure field by the density "rho" and finally by the area, to get the force in newtons (N).

On the other hand....
In the case of multiphase solvers like the ones you are using (interFoam, rasInterFoam, etc...), in the transportProperties file that is present in the "constant" directory of a case, you need to specify the density of each medium present.

In this case, the pressure value specified in the file "pd", is directly in the usual physical units of pressure (Pa = kg/ms²).

So... in the case which you tried (where the structure is completely submerged under water), I assume you used a "rho" value of 1000 in the "calcPressureForcesDict" file. This is why you got a force of 95818 N, instead of 95 N. Since the pressure field is already in "Pa", you dont need to multiply the pressure by density, which means, you can use a "rho" value of "1" instead of "1000" in "calcPressureForcesDict".


8. Now for modifying the utility to handle multiple phases..... You need to make no modifications :-)! Since the pressure field calculated in "pd" is directly in the usual physical units of pressure, you dont need to worry about the "gamma" field, etc... all you have to do, is to always set the density (rho) value in the "calcPressureForcesDict" file to unity (1.0).

9. The type of integration algorithm used.... I do not use any special integration algorithm (other than saying: integration => sum of many values :-)!). Before you started with the fluid solver, you discretised your domain by meshing it. So you basically broke down your surfaces into a large number of small faces, which are part of small cells. The solver calculates the pressure in each of these cells, and stores the value at the center of the cell.

At a boundary, when you ask for the pressure of a face, it basically gives you the value of pressure calculated by the solver for the cell to which that particular face belongs. Now... usually a patch consists of many faces. So to calculate the total force on a patch due to the pressure distribution on that patch, you get:

F_patch = Sigma(P_faces * A_faces)

Since you know the surface area of each of those tiny faces, and the pressure at each of those tiny faces, you just need to calculate the force at each of those tiny faces, and add them all up over the patch :-)!

9. As for the difference between the two equations.....

mesh.Sf() gives you the normal vector of the face in question, but, with the magnitude of that vector being equal to the surface area of that face.

Since force is a vector quantity, you need to not only calculate F = P*A, but you also need the direction in which this force is acting. In the case of the force acting on a surface due to a pressure, the direction of the force is always normal to the surface. Hence,

patchPressureForce = gSum(p.boundaryField()[forcePatchID]
* mesh.Sf().boundaryField()[forcePatchID]
);

is like hitting two birds with one stone :-) You get the force, and the direction of that force on a patch.

mesh.magSf() on the other hand, gives you only the surface area of the face in question as a scalar (basically, the magnitude of the normal vector returned by mesh.Sf()..).

I wanted to calculate the mean pressure acting on each of my patches.... since the mesh is not of a uniform size at all points, I cannot add up the pressure on each face in the patch and divide it by the number of faces of that patch. Instead, I need to calculate the area weighted average... which is basically...:

P_mean = Sigma(P_faces * A_faces) / Sigma(A_faces)

Whew... that was a long post :-)! Hope it clears up all your doubts :-)!

I wish you a nice weekend!

Philippose

[paka]

Hi Philippose!

Many very thanks for your detailed explanations and time! They are great. I promise I will reply to them with the beginning of coming week.

Meanwhile, during the weekend time, I will ask one more question.

Let me summarize:
a) The calcPressureForce tool provides the forces in 3 different directions.

b) The tool does not include transient pressures yet.

In my case, the flow is horizontal - in 1-direction, in a 2-D rectangular channel domain. The vertical direction is a 2nd-direction.

When I use calcPressureForce tool, I would expect two force components, in 1st and 2nd direction. However, the horizontal components are always "0" (zero) in every point of analyzed rigid structure, whereas vertical components has values as expected.

Is it right to have horizontal components equal to zero? Intuition tells me there should be some values, at least small.

Structure has "wall" definitions at its boundaries.

Regards and thanks!
Krystian

[philippose]

Hi Krystian!

I have a few minutes (waiting for a mesh to finish) before jumping into bed... so shall try to answer you questions in the meantime :-)!


a) Yes... calcPressureForces gives you the forces in 3 directions, which are by default (as setup in the dictionary I sent you):

(1 0 0) => x-direction
(0 1 0) => y-direction
(0 0 1) => z-direction (calculated internally)

You can always change the first two vectors if you need to align the vectors along which the forces are resolved, to your system of interest.


b) I do not understand what you mean by whether the "tool includes" transient pressures or not. Basically, when you do a "transient" simulation, you save a number of time-steps during the simulation, which (if fine enough), will record your transient pressure changes.

Now... if you run calcPressureForces after such a transient simulation, it will calculate the forces on the patches for each time step, hence capturing any changes in the forces with time.

So... what exactly do you mean by support for "transient pressures" ?


c) When you say that your flow is "horizontal", do you mean that the flow is in the global "x-direction" ? Since calcPressureForces calculates the forces in the global co-ordinate system, the force components you see for a given patch would depend on the orientation of the patch surfaces with respect to the global co-ordinate system, and the vectors you have provided for "direction 1" and "direction 2".

Would it be possible for you to post an image of your simulation? so that I can understand more about the system itself? It would make things a lot easier :-)!

Have a nice day!

Philippose

[paka]

Evening, your morning, Philippose

b) Sorry for inaccurate explanation. Thoughts circulate between different topics.

Actually, I meant SHEAR FORCES.

BTW. Tomorrow morning I'm going to test the force calculation on other (vertical) wall patches.

c)
Yes it is in x-direction and patch on which I calculate those forces is also in x direction.

Since, my case is 2D and patch is along x direction, the only force component I obtain is in y direction, which is force normal to the surface, and which actually makes sense with your explanation above.

Does it mean, that the only forces which are actually calculate, in case of orthogonal structure alignment, would be normal to the patch surface?

Does it also mean, that the only tangent forces to the patch (if it is orthogonal) would come from the shear components?

Talking about a screenshot, please forgive me, I can provide you a glance of my structure, however here, we've still got Sunday so I do not have access to my work files, so if something still won't be clear till tomorrow, I promise I will include one of the screenshots tomorrow.

d) So many points
Would it be hard to add shear forces computation to the code?

I will get back to you tomorrow, Monday.

Many thanks,
Krystian

[philippose]

Hi Krystian!

A Good day to you :-)!

No problem about the confusion regarding "transient pressure" and "shear forces" :-)!

So, in answer to your question c):

The direction of the force acting on a surface due to pressure on that surface, is always normal to the surface. But then... you need to be a little careful here....

As you know, you are discretising your domain by means of a mesh. This means, that the patch that you are trying the calculate the force on, is actually composed of a large number of small faces or surfaces...... As I had mentioned earlier, the solver calculates the pressure in the cells to which each of those faces belong.

Now... when you calculate the force on the patch, what you are actually doing, is finding the force on each of those small faces, due to the pressure acting on that face. Each of these forces calculated, will act on the corresponding face, in a direction normal to that face. Then, you finally add up all these force vectors over all the faces on the patch, to give one final vector whose magnitude and direction is the resultant force on that patch.

This final resultant force, and its direction, depend on the shape of the patch, and the pressure distribution over that patch. So, if your patch is a flat surface parallel to the global x-axis, the only force you will see due to a pressure acting on that surface, will be in the y-axis, because all your surface normals will be parallel to the y-axis.

As for question d)....

It is not difficult to add code to calculate the shear forces on a patch. The difficulty arises, in determining what equations to use for calculating the shear forces.

When you have purely laminar flow over a surface, then the equation for shear force is the basic equation:

F_shear = Sigma(A_face * mu * (du/dh)_face)

Where mu is the dynamic viscosity, and (du/dh)_face is the gradient of the velocity at each face along a direction h, which is normal to the face.

On the other hand, when you have turbulent flow, if I am not mistaken, you also need to account for the variation in local viscosity of the medium near the walls, so you have an additional equation for the turbulent viscous forces.

You can have a look at the file "liftDrag.C", which is present in "OpenFOAM-<version>/src/postProcessing/liftDrag", to gain a idea of how these forces can be calculated.

Hope this helps :-)!

Have a nice day!

Philippose

[paka]

Philippose,

Once again many very thanks for your detailed explanations. They are loud and clear No further questions for now.

Thank again and have a wonderful day
Krystian

[mgz1985]

hi phillipose,

I am an year late to this discussion but it is very important and interesting for me. I am solving an incompresible flow for a 2d airfoil. everything is OK but i do not have the lift drag utility in my version so i was thinking to add your patch to my utilities. i believe i can make a change for incorporating the shear forces in it.

but can u tell me how to add this to a disk image of OF on a macintosh .

thanx a lot

[mbots]

Hello,

I am trying to implement the calcPressureForces utility. Should it be working for OpenFOAM version 1606+ as well?

I am running into an error while I am compiling the new function. Does someone know what I did wrong? The steps I took are:

1. copy the calcPressureForces folder to: OpenFOAM/OpenFOAM-v3.0+/src/postProcessing/functionObjects/field/

2. copy the files calcPressureForces.C and createP.H to OpenFOAM/OpenFOAM-v3.0+/src/postProcessing/functionObjects/field/lnInclude/

3. add the code from the 'files' file from CalcPressureForces/Make to the files file in field/Make.

4. execute wmake in the directory OpenFOAM-v3.0+/src/postProcessing/functionObjects/field/

5. then I run into the following error:
calcPressureForces/calcPressureForces.C: In function ‘int main(int, char**)’:
calcPressureForces/calcPressureForces.C:146:31: error: ‘endTime’ was not declared in this scope
for (label i=startTime; i<endTime; i++)


Thank you for your help,
Michiel