1 package net.sf.openrocket.aerodynamics.barrowman;
3 import static net.sf.openrocket.aerodynamics.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.WarningSet;
9 import net.sf.openrocket.rocketcomponent.BodyTube;
10 import net.sf.openrocket.rocketcomponent.RocketComponent;
11 import net.sf.openrocket.rocketcomponent.SymmetricComponent;
12 import net.sf.openrocket.rocketcomponent.Transition;
13 import net.sf.openrocket.util.Coordinate;
14 import net.sf.openrocket.util.LinearInterpolator;
15 import net.sf.openrocket.util.MathUtil;
16 import net.sf.openrocket.util.PolyInterpolator;
21 * Calculates the aerodynamic properties of a <code>SymmetricComponent</code>.
23 * CP and CNa are calculated by the Barrowman method extended to account for body lift
24 * by the method presented by Galejs. Supersonic CNa and CP are assumed to be the
25 * same as the subsonic values.
28 * @author Sampo Niskanen <sampo.niskanen@iki.fi>
30 public class SymmetricComponentCalc extends RocketComponentCalc {
32 public static final double BODY_LIFT_K = 1.1;
34 private final SymmetricComponent component;
36 private final double length;
37 private final double r1, r2;
38 private final double fineness;
39 private final Transition.Shape shape;
40 private final double param;
41 private final double area;
43 public SymmetricComponentCalc(RocketComponent c) {
45 if (!(c instanceof SymmetricComponent)) {
46 throw new IllegalArgumentException("Illegal component type "+c);
48 this.component = (SymmetricComponent) c;
51 length = component.getLength();
52 r1 = component.getForeRadius();
53 r2 = component.getAftRadius();
55 fineness = length / (2*Math.abs(r2-r1));
57 if (component instanceof BodyTube) {
61 } else if (component instanceof Transition) {
62 shape = ((Transition)component).getType();
63 param = ((Transition)component).getShapeParameter();
64 area = Math.abs(Math.PI * (r1*r1 - r2*r2));
66 throw new UnsupportedOperationException("Unknown component type " +
67 component.getComponentName());
72 private boolean isTube = false;
73 private double cnaCache = Double.NaN;
74 private double cpCache = Double.NaN;
78 * Calculates the non-axial forces produced by the fins (normal and side forces,
79 * pitch, yaw and roll moments, CP position, CNa).
81 * This method uses the Barrowman method for CP and CNa calculation and the
82 * extension presented by Galejs for the effect of body lift.
84 * The CP and CNa at supersonic speeds are assumed to be the same as those at
88 public void calculateNonaxialForces(FlightConditions conditions,
89 AerodynamicForces forces, WarningSet warnings) {
91 // Pre-calculate and store the results
92 if (Double.isNaN(cnaCache)) {
93 final double r0 = component.getForeRadius();
94 final double r1 = component.getAftRadius();
96 if (MathUtil.equals(r0, r1)) {
102 final double A0 = Math.PI * pow2(r0);
103 final double A1 = Math.PI * pow2(r1);
105 cnaCache = 2 * (A1 - A0);
106 System.out.println("cnaCache = "+cnaCache);
107 cpCache = (component.getLength() * A1 - component.getFullVolume()) / (A1 - A0);
113 // If fore == aft, only body lift is encountered
115 cp = getLiftCP(conditions, warnings);
117 cp = new Coordinate(cpCache,0,0,cnaCache * conditions.getSincAOA() /
118 conditions.getRefArea()).average(getLiftCP(conditions,warnings));
122 forces.CNa = cp.weight;
123 forces.CN = forces.CNa * conditions.getAOA();
124 forces.Cm = forces.CN * cp.x / conditions.getRefLength();
126 forces.CrollDamp = 0;
127 forces.CrollForce = 0;
132 // Add warning on supersonic flight
133 if (conditions.getMach() > 1.1) {
134 warnings.add("Body calculations may not be entirely accurate at supersonic speeds.");
142 * Calculate the body lift effect according to Galejs.
144 protected Coordinate getLiftCP(FlightConditions conditions, WarningSet warnings) {
145 double area = component.getComponentPlanformArea();
146 double center = component.getComponentPlanformCenter();
149 * Without this extra multiplier the rocket may become unstable at apogee
150 * when turning around, and begin oscillating horizontally. During the flight
151 * of the rocket this has no effect. It is effective only when AOA > 45 deg
152 * and the velocity is less than 15 m/s.
155 if ((conditions.getMach() < 0.05) && (conditions.getAOA() > Math.PI/4)) {
156 mul = pow2(conditions.getMach() / 0.05);
159 return new Coordinate(center, 0, 0, mul*BODY_LIFT_K * area/conditions.getRefArea() *
160 conditions.getSinAOA() * conditions.getSincAOA()); // sin(aoa)^2 / aoa
165 private LinearInterpolator interpolator = null;
168 public double calculatePressureDragForce(FlightConditions conditions,
169 double stagnationCD, double baseCD, WarningSet warnings) {
171 if (component instanceof BodyTube)
174 if (!(component instanceof Transition)) {
175 throw new RuntimeException("Pressure calculation of unknown type: "+
176 component.getComponentName());
179 // Check for simple cases first
183 if (length < 0.001) {
185 return stagnationCD * area / conditions.getRefArea();
187 return baseCD * area / conditions.getRefArea();
192 // Boattail drag computed directly from base drag
196 double cd = baseCD * area / conditions.getRefArea();
199 return cd * (3-fineness)/2;
203 assert(r1 < r2); // Tube and boattail have been checked already
206 // All nose cones and shoulders from pre-calculated and interpolating
207 if (interpolator == null) {
208 calculateNoseInterpolator();
211 return interpolator.getValue(conditions.getMach()) * area / conditions.getRefArea();
217 * Experimental values of pressure drag for different nose cone shapes with a fineness
218 * ratio of 3. The data is taken from 'Collection of Zero-Lift Drag Data on Bodies
219 * of Revolution from Free-Flight Investigations', NASA TR-R-100, NTRS 19630004995,
222 * This data is extrapolated for other fineness ratios.
225 private static final LinearInterpolator ellipsoidInterpolator = new LinearInterpolator(
226 new double[] { 1.2, 1.25, 1.3, 1.4, 1.6, 2.0, 2.4 },
227 new double[] {0.110, 0.128, 0.140, 0.148, 0.152, 0.159, 0.162 /* constant */ }
229 private static final LinearInterpolator x14Interpolator = new LinearInterpolator(
230 new double[] { 1.2, 1.3, 1.4, 1.6, 1.8, 2.2, 2.6, 3.0, 3.6},
231 new double[] {0.140, 0.156, 0.169, 0.192, 0.206, 0.227, 0.241, 0.249, 0.252}
233 private static final LinearInterpolator x12Interpolator = new LinearInterpolator(
234 new double[] {0.925, 0.95, 1.0, 1.05, 1.1, 1.2, 1.3, 1.7, 2.0},
235 new double[] { 0, 0.014, 0.050, 0.060, 0.059, 0.081, 0.084, 0.085, 0.078}
237 private static final LinearInterpolator x34Interpolator = new LinearInterpolator(
238 new double[] { 0.8, 0.9, 1.0, 1.06, 1.2, 1.4, 1.6, 2.0, 2.8, 3.4},
239 new double[] { 0, 0.015, 0.078, 0.121, 0.110, 0.098, 0.090, 0.084, 0.078, 0.074}
241 private static final LinearInterpolator vonKarmanInterpolator = new LinearInterpolator(
242 new double[] { 0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4, 1.6, 2.0, 3.0},
243 new double[] { 0, 0.010, 0.027, 0.055, 0.070, 0.081, 0.095, 0.097, 0.091, 0.083}
245 private static final LinearInterpolator lvHaackInterpolator = new LinearInterpolator(
246 new double[] { 0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4, 1.6, 2.0 },
247 new double[] { 0, 0.010, 0.024, 0.066, 0.084, 0.100, 0.114, 0.117, 0.113 }
249 private static final LinearInterpolator parabolicInterpolator = new LinearInterpolator(
250 new double[] {0.95, 0.975, 1.0, 1.05, 1.1, 1.2, 1.4, 1.7},
251 new double[] { 0, 0.016, 0.041, 0.092, 0.109, 0.119, 0.113, 0.108}
253 private static final LinearInterpolator parabolic12Interpolator = new LinearInterpolator(
254 new double[] { 0.8, 0.9, 0.95, 1.0, 1.05, 1.1, 1.3, 1.5, 1.8},
255 new double[] { 0, 0.016, 0.042, 0.100, 0.126, 0.125, 0.100, 0.090, 0.088}
257 private static final LinearInterpolator parabolic34Interpolator = new LinearInterpolator(
258 new double[] { 0.9, 0.95, 1.0, 1.05, 1.1, 1.2, 1.4, 1.7},
259 new double[] { 0, 0.023, 0.073, 0.098, 0.107, 0.106, 0.089, 0.082}
261 private static final LinearInterpolator bluntInterpolator = new LinearInterpolator();
263 for (double m=0; m<3; m+=0.05)
264 bluntInterpolator.addPoint(m, BarrowmanCalculator.calculateStagnationCD(m));
268 * Calculate the LinearInterpolator 'interpolator'. After this call, if can be used
269 * to get the pressure drag coefficient at any Mach number.
271 * First, the transonic/supersonic region is computed. For conical and ogive shapes
272 * this is calculated directly. For other shapes, the values for fineness-ratio 3
273 * transitions are taken from the experimental values stored above (for parameterized
274 * shapes the values are interpolated between the parameter values). These are then
275 * extrapolated to the current fineness ratio.
277 * Finally, if the first data points in the interpolator are not zero, the subsonic
278 * region is interpolated in the form Cd = a*M^b + Cd(M=0).
280 @SuppressWarnings("null")
281 private void calculateNoseInterpolator() {
282 LinearInterpolator int1=null, int2=null;
285 interpolator = new LinearInterpolator();
287 double r = component.getRadius(0.99*length);
288 double sinphi = (r2-r)/MathUtil.hypot(r2-r, 0.01*length);
291 * Take into account nose cone shape. Conical and ogive generate the interpolator
292 * directly. Others store a interpolator for fineness ratio 3 into int1, or
293 * for parameterized shapes store the bounding fineness ratio 3 interpolators into
294 * int1 and int2 and set 0 <= p <= 1 according to the bounds.
298 interpolator = calculateOgiveNoseInterpolator(0, sinphi); // param==0 -> conical
302 interpolator = calculateOgiveNoseInterpolator(param, sinphi);
306 int1 = ellipsoidInterpolator;
311 int1 = bluntInterpolator;
312 int2 = x14Interpolator;
314 } else if (param <= 0.5) {
315 int1 = x14Interpolator;
316 int2 = x12Interpolator;
318 } else if (param <= 0.75) {
319 int1 = x12Interpolator;
320 int2 = x34Interpolator;
323 int1 = x34Interpolator;
324 int2 = calculateOgiveNoseInterpolator(0, 1/Math.sqrt(1+4*pow2(fineness)));
331 int1 = calculateOgiveNoseInterpolator(0, 1/Math.sqrt(1+4*pow2(fineness)));
332 int2 = parabolic12Interpolator;
334 } else if (param <= 0.75) {
335 int1 = parabolic12Interpolator;
336 int2 = parabolic34Interpolator;
339 int1 = parabolic34Interpolator;
340 int2 = parabolicInterpolator;
346 int1 = vonKarmanInterpolator;
347 int2 = lvHaackInterpolator;
352 throw new UnsupportedOperationException("Unknown transition shape: "+shape);
359 // Check for parameterized shape and interpolate if necessary
361 LinearInterpolator int3 = new LinearInterpolator();
362 for (double m: int1.getXPoints()) {
363 int3.addPoint(m, p*int2.getValue(m) + (1-p)*int1.getValue(m));
365 for (double m: int2.getXPoints()) {
366 int3.addPoint(m, p*int2.getValue(m) + (1-p)*int1.getValue(m));
371 // Extrapolate for fineness ratio if necessary
373 double log4 = Math.log(fineness+1) / Math.log(4);
374 for (double m: int1.getXPoints()) {
375 double stag = bluntInterpolator.getValue(m);
376 interpolator.addPoint(m, stag*Math.pow(int1.getValue(m)/stag, log4));
382 * Now the transonic/supersonic region is ok. We still need to interpolate
383 * the subsonic region, if the values are non-zero.
386 double min = interpolator.getXPoints()[0];
387 double minValue = interpolator.getValue(min);
388 if (minValue < 0.001) {
389 // No interpolation necessary
393 double cdMach0 = 0.8 * pow2(sinphi);
394 double minDeriv = (interpolator.getValue(min+0.01) - minValue)/0.01;
396 // These should not occur, but might cause havoc for the interpolation
397 if ((cdMach0 >= minValue-0.01) || (minDeriv <= 0.01)) {
401 // Cd = a*M^b + cdMach0
402 double a = minValue - cdMach0;
403 double b = minDeriv / a;
405 for (double m=0; m < minValue; m+= 0.05) {
406 interpolator.addPoint(m, a*Math.pow(m, b) + cdMach0);
411 private static final PolyInterpolator conicalPolyInterpolator =
412 new PolyInterpolator(new double[] {1.0, 1.3}, new double[] {1.0, 1.3});
414 private static LinearInterpolator calculateOgiveNoseInterpolator(double param,
416 LinearInterpolator interpolator = new LinearInterpolator();
418 // In the range M = 1 ... 1.3 use polynomial approximation
419 double cdMach1 = 2.1*pow2(sinphi) + 0.6019*sinphi;
421 double[] poly = conicalPolyInterpolator.interpolator(
423 4/(GAMMA+1) * (1 - 0.5*cdMach1), -1.1341*sinphi
426 // Shape parameter multiplier
427 double mul = 0.72 * pow2(param-0.5) + 0.82;
429 for (double m = 1; m < 1.3001; m += 0.02) {
430 interpolator.addPoint(m, mul * PolyInterpolator.eval(m, poly));
433 // Above M = 1.3 use direct formula
434 for (double m = 1.32; m < 4; m += 0.02) {
435 interpolator.addPoint(m, mul * (2.1*pow2(sinphi) + 0.5*sinphi/Math.sqrt(m*m - 1)));