1 package net.sf.openrocket.aerodynamics.barrowman;
3 import static net.sf.openrocket.models.atmosphere.AtmosphericConditions.GAMMA;
4 import static net.sf.openrocket.util.MathUtil.pow2;
5 import net.sf.openrocket.aerodynamics.AerodynamicForces;
6 import net.sf.openrocket.aerodynamics.BarrowmanCalculator;
7 import net.sf.openrocket.aerodynamics.FlightConditions;
8 import net.sf.openrocket.aerodynamics.Warning;
9 import net.sf.openrocket.aerodynamics.WarningSet;
10 import net.sf.openrocket.rocketcomponent.BodyTube;
11 import net.sf.openrocket.rocketcomponent.RocketComponent;
12 import net.sf.openrocket.rocketcomponent.SymmetricComponent;
13 import net.sf.openrocket.rocketcomponent.Transition;
14 import net.sf.openrocket.util.BugException;
15 import net.sf.openrocket.util.Coordinate;
16 import net.sf.openrocket.util.LinearInterpolator;
17 import net.sf.openrocket.util.MathUtil;
18 import net.sf.openrocket.util.PolyInterpolator;
23 * Calculates the aerodynamic properties of a <code>SymmetricComponent</code>.
25 * CP and CNa are calculated by the Barrowman method extended to account for body lift
26 * by the method presented by Galejs. Supersonic CNa and CP are assumed to be the
27 * same as the subsonic values.
30 * @author Sampo Niskanen <sampo.niskanen@iki.fi>
32 public class SymmetricComponentCalc extends RocketComponentCalc {
34 public static final double BODY_LIFT_K = 1.1;
36 private final double length;
37 private final double foreRadius, aftRadius;
38 private final double fineness;
39 private final Transition.Shape shape;
40 private final double param;
41 private final double frontalArea;
42 private final double fullVolume;
43 private final double planformArea, planformCenter;
44 private final double sinphi;
46 public SymmetricComponentCalc(RocketComponent c) {
48 if (!(c instanceof SymmetricComponent)) {
49 throw new IllegalArgumentException("Illegal component type " + c);
51 SymmetricComponent component = (SymmetricComponent) c;
54 length = component.getLength();
55 foreRadius = component.getForeRadius();
56 aftRadius = component.getAftRadius();
58 fineness = length / (2 * Math.abs(aftRadius - foreRadius));
59 fullVolume = component.getFullVolume();
60 planformArea = component.getComponentPlanformArea();
61 planformCenter = component.getComponentPlanformCenter();
63 if (component instanceof BodyTube) {
68 } else if (component instanceof Transition) {
69 shape = ((Transition) component).getType();
70 param = ((Transition) component).getShapeParameter();
71 frontalArea = Math.abs(Math.PI * (foreRadius * foreRadius - aftRadius * aftRadius));
73 double r = component.getRadius(0.99 * length);
74 sinphi = (aftRadius - r) / MathUtil.hypot(aftRadius - r, 0.01 * length);
76 throw new UnsupportedOperationException("Unknown component type " +
77 component.getComponentName());
82 private boolean isTube = false;
83 private double cnaCache = Double.NaN;
84 private double cpCache = Double.NaN;
88 * Calculates the non-axial forces produced by the fins (normal and side forces,
89 * pitch, yaw and roll moments, CP position, CNa).
91 * This method uses the Barrowman method for CP and CNa calculation and the
92 * extension presented by Galejs for the effect of body lift.
94 * The CP and CNa at supersonic speeds are assumed to be the same as those at
98 public void calculateNonaxialForces(FlightConditions conditions,
99 AerodynamicForces forces, WarningSet warnings) {
101 // Pre-calculate and store the results
102 if (Double.isNaN(cnaCache)) {
103 final double r0 = foreRadius;
104 final double r1 = aftRadius;
106 if (MathUtil.equals(r0, r1)) {
112 final double A0 = Math.PI * pow2(r0);
113 final double A1 = Math.PI * pow2(r1);
115 cnaCache = 2 * (A1 - A0);
116 // System.out.println("cnaCache = " + cnaCache);
117 cpCache = (length * A1 - fullVolume) / (A1 - A0);
123 // If fore == aft, only body lift is encountered
125 cp = getLiftCP(conditions, warnings);
127 cp = new Coordinate(cpCache, 0, 0, cnaCache * conditions.getSincAOA() /
128 conditions.getRefArea()).average(getLiftCP(conditions, warnings));
132 forces.setCNa(cp.weight);
133 forces.setCN(forces.getCNa() * conditions.getAOA());
134 forces.setCm(forces.getCN() * cp.x / conditions.getRefLength());
136 forces.setCrollDamp(0);
137 forces.setCrollForce(0);
142 // Add warning on supersonic flight
143 if (conditions.getMach() > 1.1) {
144 warnings.add(Warning.SUPERSONIC);
152 * Calculate the body lift effect according to Galejs.
154 protected Coordinate getLiftCP(FlightConditions conditions, WarningSet warnings) {
157 * Without this extra multiplier the rocket may become unstable at apogee
158 * when turning around, and begin oscillating horizontally. During the flight
159 * of the rocket this has no effect. It is effective only when AOA > 45 deg
160 * and the velocity is less than 15 m/s.
162 * TODO: MEDIUM: This causes an anomaly to the flight results with the CP jumping at apogee
165 if ((conditions.getMach() < 0.05) && (conditions.getAOA() > Math.PI / 4)) {
166 mul = pow2(conditions.getMach() / 0.05);
169 return new Coordinate(planformCenter, 0, 0, mul * BODY_LIFT_K * planformArea / conditions.getRefArea() *
170 conditions.getSinAOA() * conditions.getSincAOA()); // sin(aoa)^2 / aoa
175 private LinearInterpolator interpolator = null;
178 public double calculatePressureDragForce(FlightConditions conditions,
179 double stagnationCD, double baseCD, WarningSet warnings) {
181 // Check for simple cases first
182 if (MathUtil.equals(foreRadius, aftRadius))
185 if (length < 0.001) {
186 if (foreRadius < aftRadius) {
187 return stagnationCD * frontalArea / conditions.getRefArea();
189 return baseCD * frontalArea / conditions.getRefArea();
194 // Boattail drag computed directly from base drag
195 if (aftRadius < foreRadius) {
198 double cd = baseCD * frontalArea / conditions.getRefArea();
201 return cd * (3 - fineness) / 2;
205 // All nose cones and shoulders from pre-calculated and interpolating
206 if (interpolator == null) {
207 calculateNoseInterpolator();
210 return interpolator.getValue(conditions.getMach()) * frontalArea / conditions.getRefArea();
216 * Experimental values of pressure drag for different nose cone shapes with a fineness
217 * ratio of 3. The data is taken from 'Collection of Zero-Lift Drag Data on Bodies
218 * of Revolution from Free-Flight Investigations', NASA TR-R-100, NTRS 19630004995,
221 * This data is extrapolated for other fineness ratios.
224 private static final LinearInterpolator ellipsoidInterpolator = new LinearInterpolator(
225 new double[] { 1.2, 1.25, 1.3, 1.4, 1.6, 2.0, 2.4 },
226 new double[] { 0.110, 0.128, 0.140, 0.148, 0.152, 0.159, 0.162 /* constant */}
228 private static final LinearInterpolator x14Interpolator = new LinearInterpolator(
229 new double[] { 1.2, 1.3, 1.4, 1.6, 1.8, 2.2, 2.6, 3.0, 3.6 },
230 new double[] { 0.140, 0.156, 0.169, 0.192, 0.206, 0.227, 0.241, 0.249, 0.252 }
232 private static final LinearInterpolator x12Interpolator = new LinearInterpolator(
233 new double[] { 0.925, 0.95, 1.0, 1.05, 1.1, 1.2, 1.3, 1.7, 2.0 },
234 new double[] { 0, 0.014, 0.050, 0.060, 0.059, 0.081, 0.084, 0.085, 0.078 }
236 private static final LinearInterpolator x34Interpolator = new LinearInterpolator(
237 new double[] { 0.8, 0.9, 1.0, 1.06, 1.2, 1.4, 1.6, 2.0, 2.8, 3.4 },
238 new double[] { 0, 0.015, 0.078, 0.121, 0.110, 0.098, 0.090, 0.084, 0.078, 0.074 }
240 private static final LinearInterpolator vonKarmanInterpolator = new LinearInterpolator(
241 new double[] { 0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4, 1.6, 2.0, 3.0 },
242 new double[] { 0, 0.010, 0.027, 0.055, 0.070, 0.081, 0.095, 0.097, 0.091, 0.083 }
244 private static final LinearInterpolator lvHaackInterpolator = new LinearInterpolator(
245 new double[] { 0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4, 1.6, 2.0 },
246 new double[] { 0, 0.010, 0.024, 0.066, 0.084, 0.100, 0.114, 0.117, 0.113 }
248 private static final LinearInterpolator parabolicInterpolator = new LinearInterpolator(
249 new double[] { 0.95, 0.975, 1.0, 1.05, 1.1, 1.2, 1.4, 1.7 },
250 new double[] { 0, 0.016, 0.041, 0.092, 0.109, 0.119, 0.113, 0.108 }
252 private static final LinearInterpolator parabolic12Interpolator = new LinearInterpolator(
253 new double[] { 0.8, 0.9, 0.95, 1.0, 1.05, 1.1, 1.3, 1.5, 1.8 },
254 new double[] { 0, 0.016, 0.042, 0.100, 0.126, 0.125, 0.100, 0.090, 0.088 }
256 private static final LinearInterpolator parabolic34Interpolator = new LinearInterpolator(
257 new double[] { 0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4, 1.7 },
258 new double[] { 0, 0.023, 0.073, 0.098, 0.107, 0.106, 0.089, 0.082 }
260 private static final LinearInterpolator bluntInterpolator = new LinearInterpolator();
262 for (double m = 0; m < 3; m += 0.05)
263 bluntInterpolator.addPoint(m, BarrowmanCalculator.calculateStagnationCD(m));
267 * Calculate the LinearInterpolator 'interpolator'. After this call, if can be used
268 * to get the pressure drag coefficient at any Mach number.
270 * First, the transonic/supersonic region is computed. For conical and ogive shapes
271 * this is calculated directly. For other shapes, the values for fineness-ratio 3
272 * transitions are taken from the experimental values stored above (for parameterized
273 * shapes the values are interpolated between the parameter values). These are then
274 * extrapolated to the current fineness ratio.
276 * Finally, if the first data points in the interpolator are not zero, the subsonic
277 * region is interpolated in the form Cd = a*M^b + Cd(M=0).
279 @SuppressWarnings("null")
280 private void calculateNoseInterpolator() {
281 LinearInterpolator int1 = null, int2 = null;
284 interpolator = new LinearInterpolator();
288 * Take into account nose cone shape. Conical and ogive generate the interpolator
289 * directly. Others store a interpolator for fineness ratio 3 into int1, or
290 * for parameterized shapes store the bounding fineness ratio 3 interpolators into
291 * int1 and int2 and set 0 <= p <= 1 according to the bounds.
295 interpolator = calculateOgiveNoseInterpolator(0, sinphi); // param==0 -> conical
299 interpolator = calculateOgiveNoseInterpolator(param, sinphi);
303 int1 = ellipsoidInterpolator;
308 int1 = bluntInterpolator;
309 int2 = x14Interpolator;
311 } else if (param <= 0.5) {
312 int1 = x14Interpolator;
313 int2 = x12Interpolator;
314 p = (param - 0.25) * 4;
315 } else if (param <= 0.75) {
316 int1 = x12Interpolator;
317 int2 = x34Interpolator;
318 p = (param - 0.5) * 4;
320 int1 = x34Interpolator;
321 int2 = calculateOgiveNoseInterpolator(0, 1 / MathUtil.safeSqrt(1 + 4 * pow2(fineness)));
322 p = (param - 0.75) * 4;
328 int1 = calculateOgiveNoseInterpolator(0, 1 / MathUtil.safeSqrt(1 + 4 * pow2(fineness)));
329 int2 = parabolic12Interpolator;
331 } else if (param <= 0.75) {
332 int1 = parabolic12Interpolator;
333 int2 = parabolic34Interpolator;
334 p = (param - 0.5) * 4;
336 int1 = parabolic34Interpolator;
337 int2 = parabolicInterpolator;
338 p = (param - 0.75) * 4;
343 int1 = vonKarmanInterpolator;
344 int2 = lvHaackInterpolator;
349 throw new UnsupportedOperationException("Unknown transition shape: " + shape);
352 if (p < 0 || p > 1.00001) {
353 throw new BugException("Inconsistent parameter value p=" + p + " shape=" + shape);
357 // Check for parameterized shape and interpolate if necessary
359 LinearInterpolator int3 = new LinearInterpolator();
360 for (double m : int1.getXPoints()) {
361 int3.addPoint(m, p * int2.getValue(m) + (1 - p) * int1.getValue(m));
363 for (double m : int2.getXPoints()) {
364 int3.addPoint(m, p * int2.getValue(m) + (1 - p) * int1.getValue(m));
369 // Extrapolate for fineness ratio if necessary
371 double log4 = Math.log(fineness + 1) / Math.log(4);
372 for (double m : int1.getXPoints()) {
373 double stag = bluntInterpolator.getValue(m);
374 interpolator.addPoint(m, stag * Math.pow(int1.getValue(m) / stag, log4));
380 * Now the transonic/supersonic region is ok. We still need to interpolate
381 * the subsonic region, if the values are non-zero.
384 double min = interpolator.getXPoints()[0];
385 double minValue = interpolator.getValue(min);
386 if (minValue < 0.001) {
387 // No interpolation necessary
391 double cdMach0 = 0.8 * pow2(sinphi);
392 double minDeriv = (interpolator.getValue(min + 0.01) - minValue) / 0.01;
394 // These should not occur, but might cause havoc for the interpolation
395 if ((cdMach0 >= minValue - 0.01) || (minDeriv <= 0.01)) {
399 // Cd = a*M^b + cdMach0
400 double a = minValue - cdMach0;
401 double b = minDeriv / a;
403 for (double m = 0; m < minValue; m += 0.05) {
404 interpolator.addPoint(m, a * Math.pow(m, b) + cdMach0);
409 private static final PolyInterpolator conicalPolyInterpolator =
410 new PolyInterpolator(new double[] { 1.0, 1.3 }, new double[] { 1.0, 1.3 });
412 private static LinearInterpolator calculateOgiveNoseInterpolator(double param,
414 LinearInterpolator interpolator = new LinearInterpolator();
416 // In the range M = 1 ... 1.3 use polynomial approximation
417 double cdMach1 = 2.1 * pow2(sinphi) + 0.6019 * sinphi;
419 double[] poly = conicalPolyInterpolator.interpolator(
420 1.0 * sinphi, cdMach1,
421 4 / (GAMMA + 1) * (1 - 0.5 * cdMach1), -1.1341 * sinphi
424 // Shape parameter multiplier
425 double mul = 0.72 * pow2(param - 0.5) + 0.82;
427 for (double m = 1; m < 1.3001; m += 0.02) {
428 interpolator.addPoint(m, mul * PolyInterpolator.eval(m, poly));
431 // Above M = 1.3 use direct formula
432 for (double m = 1.32; m < 4; m += 0.02) {
433 interpolator.addPoint(m, mul * (2.1 * pow2(sinphi) + 0.5 * sinphi / MathUtil.safeSqrt(m * m - 1)));