[Doug68]

Hi All,

I'm running automotive type aero simulations and I wish to work out the center of pressure for the vehicle.
I've spent a lot of time search on how this might be done but have had no luck with this forum or anywhere else on the web that I know of.

Which means the solution is either really hard or really easy I'm prepared to believe its the second option, but all I've managed to come up with is a time consuming method for working this out.

I'll explain what I am doing here as you might spot a hole in my math or better yet come up with a less long winded way of solving the problem.

In the diagram we know My (moment about y axis), Fx and Fz (forces in x and z directions) from forces.dat.
At any angle 'A' we can then work out the force (Fm) in the moment. If you run through the possible angles though this gives you a curve and not a point.
Therefore the simulation then needs to be run again with CofR set to a point reasonably past where we would estimate CoP to be (the other side of it flow wise from the first run).
This gives us another set of data to create a curve from, the curves cross twice, but one of the locations will probably look unreasonable.
The attached curves are from doing the above in excel.

I guess the final verification is to run the simulation a third time with the CoR at the calculate location and see that the Y moment is close to zero.

[julien.decharentenay]

Hi Doug,

I went through a similar exercise recently. At the time I noticed that the forces and moment output in the forces.dat file were not orthogonal and I think that they should be.

I went around the problem by assessing the center of pressure position in ParaView. But it is limited to the pressure induced forces and does not include viscous forces.

In ParaView I did the following:
- load the result on the vehicle surface only;
- enter the following formula in a calculator:

forces = p*1.205*Normals // I used simpleFoam so the pressure is normalised
Moment_downforce = forces * Pos.y() // Chech what Pos.y() is as I am typing from memory

- Apply surface integral on the above
- The y-position of the downforce is calculated as Int(Moment_downforce)/Int(area)

The above is from memory - and I am still playing around with it, but I think that the overall idea is there.

Let me know if you see any massive hole or discover an alternative solution.

[Doug68]

Hi Julien,

Can I ask you to explain:

"I noticed that the forces and moment output in the forces.dat file were not orthogonal"

A bit further, that's a bit of a worry.
Did you mean that say the moment about Y is not taken as normal to the XY or ZY plane or something else?

Thanks.

[julien.decharentenay]

Hi Doug,

I am not sure I did not get to the bottom of it.

In regards to my comment, I will try to explain. First: [Disclaimer] geometry has never been my strong point.

My understanding is that the moment is:

Moment = Vector(CoR,CoP) x Forces

Where Vector(CoR,CoP) is the vector starting at the Center of Reference (nominated by the user) and the Center of Pressure (more correctly center of application of forces) and "x" is the cross product.

Consequently to the moment being the cross product of the forces, one should have (I believe):

Moment.Forces = 0
Where "." is the dot product.

In my simulation I had the following outputs in the forces.dat:
Forces (pressure only): (47.4134 5155.76 -8312.7)
Moments (pressure only): (-17641.3 -45.1214 -155.2)

Moments.Forces = -47.4*17641.3 - 5155.7*45.1 + 8312.7*155.2 = 221,061

Is my geometry failing me?

Kind regards,
Julien

[GPan1]

I'd also be interested if you have got this to work.

One method I tried was using these equations

Mx=-Fy*z+Fz*y
My=Fx*z-Fz*x
Mz=-Fx*y+Fy*x

However, I ended up with some ridiculous answer far far away from my geometry. So I tried exporting the pressure value at each point (also the corresponding co-ordinate), multiplied each x y z length by the pressure, got the average for the whole body and divided by the average pressure which gave me a reasonable result, however I feel I may have calculated centre of geometry instead of CP. Was my idea flawed?

[Doug68]

I've not got around to chasing through Julien's information yet.
I think the method I have gets quite close, but I am not sure it is 100% accurate.

Its kind of interesting to me that this subject is not better explained or developed within the software. I would have thought that for any sort of moving body type analysis that knowing the centre of pressure would be a common thing to want to know?

[GPan1]

Your method makes sense, however in my case I don't get complete convergence and so I have had to run forces and forceCoeffs over a range where my residuals oscillate close to a value and take an average. So running the simulation three times for each time I change values of pitch/roll/yaw on the car is not practical for me.

Unless I can run several force functions, each with a different CofR?

I agree, I'm sure commercial packages do calculate CP for you, it's strange why it's not included (or easily found) in OpenFoam, and that there's little to no documentation on it.

[julien.decharentenay]

It would be great if we could square off the issue and find some robust way to do it.

Gpan1: I am curious how you went to inverse the Moment system. Whenever I tried to apply substitution to it, I could not close the system. In the approach you propose, do you use the surface normal at any point?

Doug, do you know if your approach could be applied through swak4Foam as a post-processing step? It should be reasonably straightforward to calculate moment for multiple CoR in one go. What do you think?

[GPan1]

That method is from another thread on here:

http://www.cfd-online.com/Forums/openfoam-post-processing/87224-center-pressure-aero-body.html

I'm not sure what you mean by inverting it - from what I understood, you have F* and M* for all directions and you manipulate the equations to give an x,y,z co-ordinate (which would be the CP). I think.

[julien.decharentenay]

Thanks for the link. I was meaning that the system M = F^X is a system with 2 equations and 3 unknown as the equations are not independent from each other (see post from Greg Givogue, 9 May 2011 from the link you provided).

I have been reading the posts at http://www.cfd-online.com/Forums/mai...tml#post316275 a couple of time and I am still trying to get a good grasp on it. When I can get spare time, I will check it up using ParaView as it saves from having a run a simulation everytime. Once I get the method worked up in ParaView, I will try to translate it back to openFoam.

My understanding so far:
- the system M = F^X is not sufficient to extract the location of CoP;
- One need to renormalise the system in M - MxF F = F^X - I don't quite understand why as I though that M = F^X would imply that MxF = 0; and
- One also need to impose a additional/arbitrary constraint (0 = FxX for example) to define a plane where the CoP is located.

Let me know if you have any further success.

[Doug68]

Originally Posted by julien.decharentenay

Doug, do you know if your approach could be applied through swak4Foam as a post-processing step? It should be reasonably straightforward to calculate moment for multiple CoR in one go. What do you think?

Sorry I have no experience with swak4foam, I've had no reason to go look at it before now.

[GPan1]

@ julien - well if the model is symmetric or almost symmetric, like a car, I think you could assume the CP y position is zero (x is longitudinal direction of car), but for asymmetric body, I guess you'd have to try an iterative approach.

I found you can use Doug's method, and only run the simulation once - basically in the controlDict file where you have your functions listed, you can add multiples of the forces dictionary. i.e.


{
....
#include "forces"
#include "forces2"
...
#include "forcesn"
}


Then in your system folder, you can simple duplicate the forces dictionary, change your CofR in each one, rename it, and rename the block inside so it reads

forcesn
{
type forces
...
}

you also need to create empty folders in your case directory to have
forcesn/(first number of your run i.e. if you start running forces at 0, 0, if 150 then 150)

This will give you the moments from different axis' without the need to rerun three times.

Doug - how did you get your lines to intersect, I assume you're plotting (dz=cofr offset) against (dx+cofr offset), however each plot just shows a similar curve but do not intersect.

[Doug68]

GPan, that methods great! I've tweaked it a bit to stick it all in controlDict (see attached). I'm doing a run like this right now and it seems to work just fine.

For my problems, with the model being symmetrical the point is going to be near as makes no difference on the Y plane hence I've stuck the all the CoR's on the Y.
If its a non-symmetrical then I think it would be a question of modeling the curves as 3D ellipses. The math from my first post would give the major and minor axes for the elipse, then with 3 CoR's spaced in 3D (Hopefully) the surfaces would intersect at the CoP?

[GPan1]

That should work in theory. Although I can assume symmetry for the car I'm running on. However I have a problem with intersection in 2D as well, see attached. I'm not sure how you managed to get a plot that's the inverse of the other

[Doug68]

Get my spreadsheet from:

http://guerillamotive.com/OpenFOAM/CoP.xlsx

The bit for the curves starts at row 1014 it repeats for the 2 curves.
The raw data used for the first curve are the rows up to about 1000, the raw data for the second curve is not in the sheet, but the raw data is of no interest for this bit anyway.

Hopefully from this you can see what I did to get the plots working.

[GPan1]

Ah thanks for that! I saw the only real difference was that you had a negative my moment for your second set of curves, I couldn't get a negative moment about the y axis, however having one CofR at 0 0 0, and the other at -0.8 0 0, (and taking both viscous and pressure forces into account) the curves do cross at somewhat a reasonable value. About 2/3 of the wheelbase behind the front axle (I would've thought it might be more forward due to the high pressure regions at the front - basically an open wheeler and has a flat blocky nose), the height sounds about right though.

I'll have another go later, with a different set of data points and hopefully get the same answer.

The only other method is using those moment equations, but they're completely arbitrary and are based off of your CofR point, ie. CP could be anywhere you like.

[Doug68]

Now after looking at the results from the 3 COR controldict I'm not very happy with it.
If the math tied up then the 3 curves should cross at the same point, but they don't.
With two of the COR's being at 0,0,0 and 10,0,0 they seem determined to cross at Z 0 which I don't believe for a minute. Bugger.

[GPan1]

yeah I've had similar problems, I think the correct method is using moment equations, but I can't get that to give out any sensible answers.

[GPan1]

So those equations I posted earlier should work, it's a simple moments and force diagram. Lets say you take My and the forces applied would be at the centre of pressure, so the Total Moment My = Fx*(z_pos of CP) - Fz*(x_pos of CP) (assuming z is upwards, clockwise moment is +ve) and you can repeat for Mx and Mz, and can eliminate one of two variables in an equation involving y to get either z_pos or x_pos.

However, I noticed I had large moments about the z and x axes caused by the rotating wheels, so do you think it's reasonable to neglect the wheels when finding the centre of pressure? Bearing in mind this is an open wheeler car.

[Doug68]

Hi GPan1,

I think if you've got large moments about X or Z then the chances of the CoP being on the XZ plane are slim.
Assuming the car is symmetrical on the XZ plane then I'd set about getting those moments eliminated before trying to go further with anything else.

My car has enclosed wheels and I've not worked on having them rotate yet as part of the simulation (will do eventually). For an open wheel car though I'd think you'd have to take it into account from the get go seeing as they make up such a large % of the drag of the vehicle.

I'm going to keep banging on this, keep us updated on your progress...