Excerpt from ~/OpenFOAM/OpenFO
 

Excerpt from ~/OpenFOAM/OpenFOAM-1.3/src/postProcessing/incompressible/liftDrag/liftDrag.C :

scalar qRef = 0.5*magSqr(Uinf);
scalar Fref = qRef*Aref;

vector pressureCoeff = pressureForce/Fref;
vector viscousCoeff = viscousForce/Fref;

return
(pressureCoeff + viscousCoeff)
- flowDirection*(flowDirection & (pressureCoeff + viscousCoeff));

Can somone kindly point out what the return statement is doing in the last line. Specifically, what is being subtracted?

Thanks.

It is returning the lift part
 

It is returning the lift part of the force coefficient.

i.e.
(pressureCoeff + viscousCoeff) = total force coefficient vector

flowDirection*(flowDirection & (pressureCoeff + viscousCoeff)) = the part of the force coefficient in the direction of the flow.

Thanks a lot Eugene!
 

Thanks a lot Eugene!

For a 2D case, it appears that
 

For a 2D case, it appears that the liftDrag/OpenFOAM routines use the third number specified in the latter part of the vertices section[1] in the blockMeshDict file, as the depth (which is used in the calculation of projected area, hence the diomensionless force coefficients). Changing this number, say by an order of magnitude, changes the lift/drag coefficients etc. by the same order. Someone please correct me if I am wrong.

[1] For instance the number 0.1 in the following example:

bla bla bla...

convertToMeters 1;

vertices
(
(0 0 0)
(1.16 0 0)
(1.17 0 0)
(1.515 0 0)
(0 0.015 0)
(1.16 0.015 0)
(1.17 0.015 0)
(1.515 0.015 0)
(0 0.025 0)
(1.16 0.025 0)
(1.17 0.025 0)
(1.515 0.025 0)
(0 0.04 0)
(1.16 0.04 0)
(1.17 0.04 0)
(1.515 0.04 0)
(0 0 0.1)
(1.16 0 0.1)
(1.17 0 0.1)
(1.515 0 0.1)
(0 0.015 0.1)
(1.16 0.015 0.1)
(1.17 0.015 0.1)
(1.515 0.015 0.1)
(0 0.025 0.1)
(1.16 0.025 0.1)
(1.17 0.025 0.1)
(1.515 0.025 0.1)
(0 0.04 0.1)
(1.16 0.04 0.1)
(1.17 0.04 0.1)
(1.515 0.04 0.1)
);

Despite being a very old threa
 

Despite being a very old thread, I wish to add an important observation especially for all those trying out 2D cases and using the liftDrag utility and/or its derivatives that were posted in this forum.

As mentioned above, the distance in the third direction is used to estimate the force coefficients, in particular the force on the wall patch itself. As a result even if you specify a constant projected area (say unity) when finally calculating Cd/Cl etc., the pressure used to calculate the force will be multiplied by the area of the patch as Foam sees it.

For instance, if you consider a square obstacle with dimensions D x D x L in a channel, then for a 2D case, you might create the blockMeshDict file as shown above (i.e. L = 0.1). In this case, the liftDrag calculation will assume that the area that needs to be multiplied by the pressure to get the force is D times L. Ideally then, for a 2D case, one wants L to be unity. This needs to be reflected in the blockMeshDict.

Aside notes for the checkMesh utility maintainer:
You may want to consider adding an extra check for the 2D case. If I use a unit length in the third direction and the typical grid size in the other two directions is very small compared to unity, then checkMesh reports that 1 Mesh check (very high aspect ratio) failed due to a high aspect ratio. For a 2D case, the grid sizes in the third direction should be neglected when checkMesh estimates the aspect ratio.


[madhavan@jinshi ~]$ checkMesh . franke_et_al_re_40
/*---------------------------------------------------------------------------*\
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: 1.4.1 |
| \ / A nd | Web: http://www.openfoam.org |
| \/ M anipulation | |
\*---------------------------------------------------------------------------*/

Exec : checkMesh . franke_et_al_re_40
Date : Oct 03 2007
Time : 14:22:32
Host : jinshi
PID : 26061
Root : /home/madhavan
Case : franke_et_al_re_40
Nprocs : 1
Create time

Create polyMesh for time = constant

Time = constant

Mesh stats
points: 3064640
edges: 7656160
faces: 6121120
internal faces: 3056480
cells: 1529600
boundary patches: 5
point zones: 0
face zones: 0
cell zones: 0

Number of cells of each type:
hexahedra: 1529600
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 0

Checking topology...
Boundary definition OK.
Point usage OK.
Upper triangular ordering OK.
Topological cell zip-up check OK.
Face vertices OK.
Face-face connectivity OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
Patch Faces Points Surface
ChannelWalls 3200 6404 ok (not multiply connected)
ObstacleWalls 320 640 ok (not multiply connected)
vinlet 960 1922 ok (not multiply connected)
poutlet 960 1922 ok (not multiply connected)
frontAndBack 3059200 3064640 ok (not multiply connected)

Checking geometry...
Domain bounding box: (-0.0125 -0.015 0) (0.03750000000000001 0.015 1)
Boundary openness (3.509809764395274e-17 -2.895593055626086e-15 1.174630898377996e-17) OK.
***High aspect ratio cells found, Max aspect ratio: 32000.0000000204, number of cells 1529600
<<Writing 1529600 cells with high aspect ratio to set highAspectRatioCells
Minumum face area = 9.765624999993512e-10. Maximum face area = 3.125000000001737e-05. Face area magnitudes OK.
Min volume = 9.765624999992789e-10. Max volume = 9.765625000006704e-10. Total volume = 0.001493749999965185. Cell volumes OK.
Mesh non-orthogonality Max: 8.537736462515939e-07 average: 0
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 8.881810217858113e-12 OK.
Min/max edge length = 3.124999999998268e-05 1 OK.
All angles in faces OK.
Face flatness (1 = flat, 0 = butterfly) : average = 1 min = 0.9999999999999996
All face flatness OK.

Failed 1 mesh checks.

End

I really do not agree with thi
 

I really do not agree with this: if it were so, it would be trivial to fix the width of all OpenFOAM meshes in 2-D to unit width and be done with it.

The issue stems from the way I calculate cell volume and face area, by decomposition into triangles and pyramids. In extremely stretched cells, this will accumulate the round-off errors, sometimes in a nasty way, which will mess up the accuracy in the rest of the code. Therefore, you are much better creating the mesh such that the cells in the third direction look decent (avoiding the discretisation error in surface and volume calculation) and then using a POCKET CALCULATOR to divide the forces by the width.

Hrv

Thanks for the recommendation
 

Thanks for the recommendation Dr. Jasak. I was unaware of these other issues. So to sum up, the idea is to create a 2D geometry with a decent width in the third direction (i.e. acceptable aspect ratio). I am presently comparing the pressure and viscous drag coefficient predictions with a paper. I will report the results soon.

For 3-D surfaces, high aspect
 

For 3-D surfaces, high aspect ratio cells seem inevitable when y-plus is 1 or 2, unless the wall is discritized into an unreasonably large number of very tiny pieces. My high aspect ratio cells at the wall have low skew. This appears to me to not be a problem because the OF results match wind tunnel data. I'm new to this community, so could be missing something obvious.

Hi to all, i am calculating
 

Hi to all,

i am calculating a 3D cylinder case and try to calculate the drag coefficient. The CD/Cl tool seems to work in a 2D case, but i get a negative result for my 3D case. I sue a dynSmagotinsky at Re around 200.

If i divide Cd by my depth i am close to the result i want to have, but i get the same result for Re=200 and Re=500 what is absolutly not reasenable.

Can anybody help me, i have no idea what to do.

Fabian