CoP simulation comparison, including CFD

[Buckeye]

I put my Fat Boy CFD model to the test with a center of pressure analysis for this short, stubby design.

Since we have 3D discretized numerical solution, we can easily integrate over the rocket area and define CP per the classic NASA Beginners Guide. Since the CFD includes viscous effects, I added wall shear stress magnitude into the calculation.



I think this approach make sense, and I think I calculated correctly. Regardless, the X result (Cp) did not change very much due to the addition of wall shear in my test cases.

I compared the CFD outcome to 5 other simulators. The CP locations in the figure below are color-coded to the tool used. CFD uses numerical Navier-Stokes continuum equations governing the air flow around the entire rocket, while all the others are Barrowman-based normal force, part-wise methodology.





The CFD result for CP falls in the cluster of predictions near the leading edge of the fin root. The OR + Base Drag Hack is the outlier since it radically violates Barrowman's theory. This Fat Boy case was the prototype for Levison's hack. I don't think you are doing yourself any favors by artificially moving the CP rearwards with a giant cone stuck to the back of the rocket and giving yourself a false sense of stability.

 

[Alan15578]

I put my Fat Boy CFD model to the test with a center of pressure analysis for this short, stubby design.

Since we have 3D discretized numerical solution, we can easily integrate over the rocket area and define CP per the classic NASA Beginners Guide. Since the CFD includes viscous effects, I added wall shear stress magnitude into the calculation.

View attachment 630001

I think this approach make sense, and I think I calculated correctly. Regardless, the X result (Cp) did not change very much due to the addition of wall shear in my test cases.

I compared the CFD outcome to 5 other simulators. The CP locations in the figure below are color-coded to the tool used. CFD uses numerical Navier-Stokes continuum equations governing the air flow around the entire rocket, while all the others are Barrowman-based normal force, part-wise methodology.

View attachment 630002

View attachment 630003

The CFD result for CP falls in the cluster of predictions near the leading edge of the fin root. The OR + Base Drag Hack is the outlier since it radically violates Barrowman's theory. This Fat Boy case was the prototype for Levison's hack. I don't think you are doing yourself any favors by artificially moving the CP rearwards with a giant cone stuck to the back of the rocket and giving yourself a false sense of stability.

It would be great if you could add a wind tunnel CP result as well. Still, attaboy for discrediting the "base drag hack".

 

[Buckeye]

Now things get interesting. I stretched the Fat Boy CFD model into longer configurations called Nominal Boy (L/D = 10) and Skinny Boy (L/D = 15).



Drag results. The longer body tube naturally increases the viscous drag contribution.



CP comparisons. The predictions spread out more as the body length increases, with the CFD model showing the most conservative values. Maybe this is the difference between potential flow (Barrowman) and N-S flow (CFD) starting to show up? Long skinny rockets often need more than 1 cal of stability, especially at angle of attack. Maybe CFD is pointing in that direction?






Of course, there is a lot more investigation with CFD needed here, like grid sensitivity, turbulence models, transient, supersonic, AoA, etc. With the ease of FreeCAD and CfdOF, those studies are very much possible, especially with arbitrary rocket designs that exceed Barrowman and DATCOM restrictions.

 

[Buckeye]

It would be great if you could add a wind tunnel CP result as well. Still, attaboy for discrediting the "base drag hack".

The RASAero webpage has lots of wind tunnel/flight measurements as comparisons. They are mostly transonic and supersonic though.

 

[boatgeek]

I would be very interested to see an angle of attack analysis. IIRC, when a friend did a CFD analysis for a square rocket, the CP moved an awful lot between dead ahead and 2 degrees of AoA. I’m not nearly as worried about stability with zero AoA as I am with stability with a small AoA. Granted, negative stability at zero AoA may result in somewhat squirrelly flights so you want at least a little positive stability there.

 

[bad_idea]

I don't think you are doing yourself any favors by artificially moving the CP rearwards with a giant cone stuck to the back of the rocket and giving yourself a false sense of stability.

I had suspected this might be the case. Glad to see some hard analysis. Thank you for sharing it!

Any chance you could expand on David Carter's presentation and show how to analyze flutter?

 

[Buckeye]

I had suspected this might be the case. Glad to see some hard analysis. Thank you for sharing it!

Any chance you could expand on David Carter's presentation and show how to analyze flutter?

Aero-elasticity is not my thing. Hopefully some FEA people will chime in and expand on it.

 

[kjhambrick]

@Buckeye --

Thanks for this thread and for your other thread CFD with FreeCAD!

WoW !

I've had FreeCAD installed for a while and after reading your threads I am trying to get Rocket Workbench going on my Slackware Linux Laptop.

Thanks again !

-- kjh

 

[JoePfeiffer]

Thanks, this is extremely useful information.

 

[SiboVG]

Interesting results, thanks @Buckeye! Would love to see more comparisons in the future.

 

[Jeff Lassahn]

Interesting stuff. And the results seem pretty plausible to me.

Is it relatively easy to do the CP integral calculation over subsets of the surface? That would be interesting because Barrowman calculates terms for individual parts, and integrating over e.g. just the fins or just the nose would allow comparing term by term.

My hypothesis is that the body is contributing a significant amount to the location of CP in your longer rockets, which pulls the CP forward compared to the Barrowman estimate (which ignores body tubes entirely)

 

[JoePfeiffer]

Interesting stuff. And the results seem pretty plausible to me.

Is it relatively easy to do the CP integral calculation over subsets of the surface? That would be interesting because Barrowman calculates terms for individual parts, and integrating over e.g. just the fins or just the nose would allow comparing term by term.

My hypothesis is that the body is contributing a significant amount to the location of CP in your longer rockets, which pulls the CP forward compared to the Barrowman estimate (which ignores body tubes entirely)

One of OR's extensions to Barrowman is we do consider the effect of body lift on CP.

 

[Buckeye]

Interesting stuff. And the results seem pretty plausible to me.

Is it relatively easy to do the CP integral calculation over subsets of the surface? That would be interesting because Barrowman calculates terms for individual parts, and integrating over e.g. just the fins or just the nose would allow comparing term by term.

My hypothesis is that the body is contributing a significant amount to the location of CP in your longer rockets, which pulls the CP forward compared to the Barrowman estimate (which ignores body tubes entirely)

Yes, the integral calculation can take place on any surface. It would be easiest to divide the rocket into the separate PIDs at the CAD level before CFD calculation.

Good observation on the body tube length impact.

 

[thzero]

One of OR's extensions to Barrowman is we do consider the effect of body lift on CP.

So does Rocsim... pretty easy to create a rocket, see what the CP is, then length the tube, and oh, hey the CP changed!

 

[Buckeye]

Another way to compute the CoP coordinate (x,y,z) is a moment balance.

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

There is no unique solution to this system of equations, as CoP lies on a plane. However, we can force some constraints:

Assume the CoP lies on the rocket centerline, so y=z=0. The rocket is symmetric about the y-axis, so Fy =0

However, it just dawned on me (duh) that my 3-fin rocket example is not symmetric about the z-axis, so there is a net z-force, even at 0 AoA. Thus, the moment equations reduce to:

My=-Fz*x

The force/moment summation is easily extracted from FreeCAD CfdOF workbench.

x = 228 mm

I call this "CFD-Moment", and the previous is now called "CFD-Integral." The moment calculation is 6 mm behind the pressure calculation, but still nowhere near the overly optimistic Base Drag Hack.

The moment method will work well for angle of attack simulations.

Here is an update to the results:

 

[boatgeek]

However, it just dawned on me (duh) that my 3-fin rocket example is not symmetric about the z-axis, so there is a net z-force, even at 0 AoA.

For the last week or so, I've been working on a 4-fin model to address that issue. I'm having some issues with getting a test model meshed, so I need to do some further troubleshooting. In an ideal world, I'll run a bunch of AoAs and possibly different lengths of models (current testbed is 5:1 L).

 

[Jeff Lassahn]

Another way to compute the CoP coordinate (x,y,z) is a moment balance.

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

There is no unique solution to this system of equations, as CoP lies on a plane. However, we can force some constraints:

Assume the CoP lies on the rocket centerline, so y=z=0. The rocket is symmetric about the y-axis, so Fy =0

However, it just dawned on me (duh) that my 3-fin rocket example is not symmetric about the z-axis, so there is a net z-force, even at 0 AoA. Thus, the moment equations reduce to:

My=-Fz*x

The force/moment summation is easily extracted from FreeCAD CfdOF workbench.

x = 228 mm

I call this "CFD-Moment", and the previous is now called "CFD-Integral." The moment calculation is 6 mm behind the pressure calculation, but still nowhere near the overly optimistic Base Drag Hack.

The moment method will work well for angle of attack simulations.

Here is an update to the results:

View attachment 632517

View attachment 632519

I think the basic approach of using the torque moments is a good one, but I think there's still some things going wrong somewhere...

At zero angle of attack, a model with a rotationally symmetric body and three identical fins is still symmetric enough that the total torque should be zero. The body obviously shouldn't have torque to the side because all directions are equal to it. And the fins shouldn't because each fin should produce an identical force, but 120 degrees apart. Three copies of a vector rotated by 120 degrees each in the YZ plane should form an equilateral triangle and cancel.

So if the three finned rocket is giving you a nonzero total torque that must either be a calculation error or numerical inaccuracy somewhere.

Coming up with a rock-solid definition if "Center of Pressure" that's meaningful in a full 3D simulation is actually pretty damned hard. Even in the symmetric cases where people are doing the math in 2D it's still pretty subtle.

The thing that's called "Center of Pressure" in Barroman's thesis (and by extension in OpenRocket and related stuff) isn't actually the pressure integral from e.g. the NASA link above, it's actually an integral over (The moment in the X direction computed around the tip of the nose of (The derivative vs angle of attack of (the normalized coefficient of (dynamic pressure))))

I'm thinking in the 3D case this thing is not even a number, it's a 2-tensor sort of analogous to the inertia tensor from rotational dynamics. That is, it's a matrix where you plug in the angle of attack (as an axial vector) and you get out a torque vector, but there's no reason to expect the torque to be in the same plane as the angle of attack.

Anyway all this to say, it's worth stepping back from "Coefficient of Pressure" abstractions and thinking about "When I change the angle of attack of the rocket, what's the resulting torque trying to push it back" because when the notation gets too abstract lots of weird crap can get lost in the mix.

 

[Buckeye]

Totally agree that Center of Pressure is an abstraction, especially in 3D.

Yeah, I went back and forth on the non-zero torque. It is there in the simulation. See the convergence plot in the CFD thread. The Fz is non-zero. I assumed it was real due to geometry non-symmetry in the cartesian coordinate system of the CFD model. Fy is actually non-zero also, but it is an order of magnitude smaller than Fz. It could all be a numerical artifact, though. One skewed mesh cell can trip up the forces. Also, the residuals converged to 1.0E-03 rather than my preferred 1.0E-04, but the force convergence certainly looked steady enough to call it good.

A 4-fin CFD model will help sort that out, as will angle of attack sims.

 

[Buckeye]

Here is a quick CFD run at 4 degrees angle of attack on the FatBoy. The sim took only 1 hour+ on 4 cores of parallel processing. Nice!

To keep the rocket on the same coordinate axes for force/moment calcs, I instead added the 4 degree angle to the boundary conditions. The front and bottom of the domain are velocity inlets with x and z components, and the top and rear are outlets. I also tilted the nearest mesh refinement box by 4 degrees to better capture the direction of the wake.







Here is CoP comparison of the tools. I could not figure how to adjust AoA in RockSim static analysis.



The CoP from the pressure integration did not change at all. The CFD moment calculation gives a plausible CoP location.

A full sweep of angles is doable, but the CFD approach may have to change to handle big angles and highly-separated flow that may occur.

 

[Buckeye]

For the last week or so, I've been working on a 4-fin model to address that issue. I'm having some issues with getting a test model meshed, so I need to do some further troubleshooting. In an ideal world, I'll run a bunch of AoAs and possibly different lengths of models (current testbed is 5:1 L).

Cool!

FreeCAD/CfdOF or some other software?

 

[boatgeek]

Cool!

FreeCAD/CfdOF or some other software?

CFDOF and FreeCAD, though I’m doing the modeling in Rhino and importing via IGS. I find the FreeCAD user interface really hard to use. Right now, gmsh is crashing just as it finishes the mesh so I may switch to a different mesher.

 

[Buckeye]

CFDOF and FreeCAD, though I’m doing the modeling in Rhino and importing via IGS. I find the FreeCAD user interface really hard to use. Right now, gmsh is crashing just as it finishes the mesh so I may switch to a different mesher.

I started with Snappy Hex Mesh because that's what I used at work. However, I don't think CfdOF offers enough control parameters for it, and I couldn't get a nice mesh. I am now using cfmesh, and it looks good on my simple rocket. I didn't even bother with gmesh-tets, and I haven't tried gmesh-poly.

Yeah, FreeCAD interface is awkward. However, I didn't have a CAD program of any kind on my home computer, so I was going to start fresh no matter what. I used ANSA at the workplace, which is heavily geared to industry CAE, and pricey.

FreeCAD also has the Rocket Workbench which makes nosecones and fin cans. That's how I am building out simple 3FNC and 4FNC rockets very quickly.

 

[boatgeek]

I started with Snappy Hex Mesh because that's what I used at work. However, I don't think CfdOF offers enough control parameters for it, and I couldn't get a nice mesh. I am now using cfmesh, and it looks good on my simple rocket. I didn't even bother with gmesh-tets, and I haven't tried gmesh-poly.

Yeah, FreeCAD interface is awkward. However, I didn't have a CAD program of any kind on my home computer, so I was going to start fresh no matter what. I used ANSA at the workplace, which is heavily geared to industry CAE, and pricey.

FreeCAD also has the Rocket Workbench which makes nosecones and fin cans. That's how I am building out simple 3FNC and 4FNC rockets very quickly.

I’ll have to look into Rocket Workbench. Body tubes and elliptical noses were OK, but any kind of not-ugly fins were a giant PITA.

 

[boatgeek]

Arrgh. Whether I use gmsh or cfmesh, I'm running into the same error:



I've double checked in that directory to make sure that both files exist. Did you run into anything like this?

 

[Buckeye]

Arrgh. Whether I use gmsh or cfmesh, I'm running into the same error:

View attachment 632697

I've double checked in that directory to make sure that both files exist. Did you run into anything like this?

You are using OpenFOAM v2312. I tried that too and ran into some errors. I backed down to v2212 which is the default download in Edit > Preferences> CfdOF. No problems with that version. Make sure you run the Dependency Checker.

 

[boatgeek]

Thank you! With that and finally getting over a initial boundary condition error, I've got a simple test simulation running. This is just a box with flow in one side and out on two other sides. I'll get to the good stuff later.

 

[aerostadt]

I am trying to setup the simple backward facing step tutorial and I cannot get CfdOF to accept the boundary condition on the faces. Any suggestions?

 

[boatgeek]

I am trying to setup the simple backward facing step tutorial and I cannot get CfdOF to accept the boundary condition on the faces. Any suggestions?

I had to select the face first, then click the boundary condition button. You'll also need to set an initial condition at at least one face. I found that pressure at an outflow face worked.

 

[Buckeye]

Here starts a CFD angle of attack study on the Fatboy.

For each angle, I rotated the entire domain and applied a velocity condition normal to the inlet face. The rocket remained aligned on the principal axes in order to keep the force and moment calculations simple.



Most simulations of small angles converged to 1.0E-04 residuals easily in a handful of iterations, like the 8 degrees below. For some reason, 12 degrees proved to be an outlier, showing more erratic behavior and non-symmetry in force and moments. I don't think this is real, but a mesh effect (Edit: It could be the real, stall on the fins, as noted by @kjhambrick.) For the highest angles of 25 and 30 degrees, some periodic wake oscillations are showing up in the steady state solution. In all cases, I took a time average over the last few hundred iterations for CP calculations.


CP was calculated with:

My=-Fz*x

except for the zero degree, since there should be no moments on this 3-fin axially-symmetric model (I think). For this case, I used the integral calculation in post #1. Both the moment and integral calculations gave nearly the same CP location within 1 mm for every angle up to 25 degrees, where differences became larger. I attribute this to the oscillations forming in the flow field of pressure and shear stress, and the integral was calculated on the final time step only (there is a way to average the field, but I haven't figured that out, yet.) The moment calcs were time-averaged, so I feel better about that data.

 

[Buckeye]

Here are the CP results, expressed as the change in location as the AOA is increased. The CFD outlier at 12 degrees is apparent. The CP movement is very small, just 11 mm forward (OR) and 8 mm (CFD) at 30 degrees. OR with the Base Drag Hack moves just 4 mm, dominated by the overly-powerful cone on the back.



I looked more closely at the Galejs paper, where he accounted for body lift and all angles of attack, which Barrowman "left out." His equations agree fairly well with measurements. (Not sure if these are the same body lift equations used in OR.) He simulated the FatBoy as well, shown in Figure 4B. I replicated the plot here and overlaid the CFD results. Delta CP is now measured in calibers.

The CFD, minus the 12 degree anomaly, matches the linear trend of the Galejs equations. The y-axis scale purposely runs out to 1.0 caliber to illustrate the absurdly- small movement of CP with angle of attack. This is why short, stubby rockets are very stable with low CP-CG margins.