git.gag.com Git - debian/openrocket/blob - core/src/net/sf/openrocket/aerodynamics/barrowman/SymmetricComponentCalc.java

   1 package net.sf.openrocket.aerodynamics.barrowman;
   2
   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;
  19
  20
  21
  22 /**
  23  * Calculates the aerodynamic properties of a <code>SymmetricComponent</code>.
  24  * <p>
  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.
  28  *
  29  *
  30  * @author Sampo Niskanen <sampo.niskanen@iki.fi>
  31  */
  32 public class SymmetricComponentCalc extends RocketComponentCalc {
  33
  34         public static final double BODY_LIFT_K = 1.1;
  35
  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;
  45
  46         public SymmetricComponentCalc(RocketComponent c) {
  47                 super(c);
  48                 if (!(c instanceof SymmetricComponent)) {
  49                         throw new IllegalArgumentException("Illegal component type " + c);
  50                 }
  51                 SymmetricComponent component = (SymmetricComponent) c;
  52
  53
  54                 length = component.getLength();
  55                 foreRadius = component.getForeRadius();
  56                 aftRadius = component.getAftRadius();
  57
  58                 fineness = length / (2 * Math.abs(aftRadius - foreRadius));
  59                 fullVolume = component.getFullVolume();
  60                 planformArea = component.getComponentPlanformArea();
  61                 planformCenter = component.getComponentPlanformCenter();
  62
  63                 if (component instanceof BodyTube) {
  64                         shape = null;
  65                         param = 0;
  66                         frontalArea = 0;
  67                         sinphi = 0;
  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));
  72
  73                         double r = component.getRadius(0.99 * length);
  74                         sinphi = (aftRadius - r) / MathUtil.hypot(aftRadius - r, 0.01 * length);
  75                 } else {
  76                         throw new UnsupportedOperationException("Unknown component type " +
  77                                         component.getComponentName());
  78                 }
  79         }
  80
  81
  82         private boolean isTube = false;
  83         private double cnaCache = Double.NaN;
  84         private double cpCache = Double.NaN;
  85
  86
  87         /**
  88          * Calculates the non-axial forces produced by the fins (normal and side forces,
  89          * pitch, yaw and roll moments, CP position, CNa).
  90          * <p>
  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.
  93          * <p>
  94          * The CP and CNa at supersonic speeds are assumed to be the same as those at
  95          * subsonic speeds.
  96          */
  97         @Override
  98         public void calculateNonaxialForces(FlightConditions conditions,
  99                         AerodynamicForces forces, WarningSet warnings) {
 100
 101                 // Pre-calculate and store the results
 102                 if (Double.isNaN(cnaCache)) {
 103                         final double r0 = foreRadius;
 104                         final double r1 = aftRadius;
 105
 106                         if (MathUtil.equals(r0, r1)) {
 107                                 isTube = true;
 108                                 cnaCache = 0;
 109                         } else {
 110                                 isTube = false;
 111
 112                                 final double A0 = Math.PI * pow2(r0);
 113                                 final double A1 = Math.PI * pow2(r1);
 114
 115                                 cnaCache = 2 * (A1 - A0);
 116                                 //                              System.out.println("cnaCache = " + cnaCache);
 117                                 cpCache = (length * A1 - fullVolume) / (A1 - A0);
 118                         }
 119                 }
 120
 121                 Coordinate cp;
 122
 123                 // If fore == aft, only body lift is encountered
 124                 if (isTube) {
 125                         cp = getLiftCP(conditions, warnings);
 126                 } else {
 127                         cp = new Coordinate(cpCache, 0, 0, cnaCache * conditions.getSincAOA() /
 128                                         conditions.getRefArea()).average(getLiftCP(conditions, warnings));
 129                 }
 130
 131                 forces.setCP(cp);
 132                 forces.setCNa(cp.weight);
 133                 forces.setCN(forces.getCNa() * conditions.getAOA());
 134                 forces.setCm(forces.getCN() * cp.x / conditions.getRefLength());
 135                 forces.setCroll(0);
 136                 forces.setCrollDamp(0);
 137                 forces.setCrollForce(0);
 138                 forces.setCside(0);
 139                 forces.setCyaw(0);
 140
 141
 142                 // Add warning on supersonic flight
 143                 if (conditions.getMach() > 1.1) {
 144                         warnings.add(Warning.SUPERSONIC);
 145                 }
 146
 147         }
 148
 149
 150
 151         /**
 152          * Calculate the body lift effect according to Galejs.
 153          */
 154         protected Coordinate getLiftCP(FlightConditions conditions, WarningSet warnings) {
 155
 156                 /*
 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.
 161                  *
 162                  * TODO: MEDIUM:  This causes an anomaly to the flight results with the CP jumping at apogee
 163                  */
 164                 double mul = 1;
 165                 if ((conditions.getMach() < 0.05) && (conditions.getAOA() > Math.PI / 4)) {
 166                         mul = pow2(conditions.getMach() / 0.05);
 167                 }
 168
 169                 return new Coordinate(planformCenter, 0, 0, mul * BODY_LIFT_K * planformArea / conditions.getRefArea() *
 170                                 conditions.getSinAOA() * conditions.getSincAOA()); // sin(aoa)^2 / aoa
 171         }
 172
 173
 174
 175         private LinearInterpolator interpolator = null;
 176
 177         @Override
 178         public double calculatePressureDragForce(FlightConditions conditions,
 179                         double stagnationCD, double baseCD, WarningSet warnings) {
 180
 181                 // Check for simple cases first
 182                 if (MathUtil.equals(foreRadius, aftRadius))
 183                         return 0;
 184
 185                 if (length < 0.001) {
 186                         if (foreRadius < aftRadius) {
 187                                 return stagnationCD * frontalArea / conditions.getRefArea();
 188                         } else {
 189                                 return baseCD * frontalArea / conditions.getRefArea();
 190                         }
 191                 }
 192
 193
 194                 // Boattail drag computed directly from base drag
 195                 if (aftRadius < foreRadius) {
 196                         if (fineness >= 3)
 197                                 return 0;
 198                         double cd = baseCD * frontalArea / conditions.getRefArea();
 199                         if (fineness <= 1)
 200                                 return cd;
 201                         return cd * (3 - fineness) / 2;
 202                 }
 203
 204
 205                 // All nose cones and shoulders from pre-calculated and interpolating
 206                 if (interpolator == null) {
 207                         calculateNoseInterpolator();
 208                 }
 209
 210                 return interpolator.getValue(conditions.getMach()) * frontalArea / conditions.getRefArea();
 211         }
 212
 213
 214
 215         /*
 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,
 219          * page 16.
 220          *
 221          * This data is extrapolated for other fineness ratios.
 222          */
 223
 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 */}
 227                         );
 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 }
 231                         );
 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 }
 235                         );
 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 }
 239                         );
 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 }
 243                         );
 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 }
 247                         );
 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 }
 251                         );
 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 }
 255                         );
 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 }
 259                         );
 260         private static final LinearInterpolator bluntInterpolator = new LinearInterpolator();
 261         static {
 262                 for (double m = 0; m < 3; m += 0.05)
 263                         bluntInterpolator.addPoint(m, BarrowmanCalculator.calculateStagnationCD(m));
 264         }
 265
 266         /**
 267          * Calculate the LinearInterpolator 'interpolator'.  After this call, if can be used
 268          * to get the pressure drag coefficient at any Mach number.
 269          *
 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.
 275          *
 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).
 278          */
 279         @SuppressWarnings("null")
 280         private void calculateNoseInterpolator() {
 281                 LinearInterpolator int1 = null, int2 = null;
 282                 double p = 0;
 283
 284                 interpolator = new LinearInterpolator();
 285
 286
 287                 /*
 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.
 292                  */
 293                 switch (shape) {
 294                 case CONICAL:
 295                         interpolator = calculateOgiveNoseInterpolator(0, sinphi); // param==0 -> conical
 296                         break;
 297
 298                 case OGIVE:
 299                         interpolator = calculateOgiveNoseInterpolator(param, sinphi);
 300                         break;
 301
 302                 case ELLIPSOID:
 303                         int1 = ellipsoidInterpolator;
 304                         break;
 305
 306                 case POWER:
 307                         if (param <= 0.25) {
 308                                 int1 = bluntInterpolator;
 309                                 int2 = x14Interpolator;
 310                                 p = param * 4;
 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;
 319                         } else {
 320                                 int1 = x34Interpolator;
 321                                 int2 = calculateOgiveNoseInterpolator(0, 1 / MathUtil.safeSqrt(1 + 4 * pow2(fineness)));
 322                                 p = (param - 0.75) * 4;
 323                         }
 324                         break;
 325
 326                 case PARABOLIC:
 327                         if (param <= 0.5) {
 328                                 int1 = calculateOgiveNoseInterpolator(0, 1 / MathUtil.safeSqrt(1 + 4 * pow2(fineness)));
 329                                 int2 = parabolic12Interpolator;
 330                                 p = param * 2;
 331                         } else if (param <= 0.75) {
 332                                 int1 = parabolic12Interpolator;
 333                                 int2 = parabolic34Interpolator;
 334                                 p = (param - 0.5) * 4;
 335                         } else {
 336                                 int1 = parabolic34Interpolator;
 337                                 int2 = parabolicInterpolator;
 338                                 p = (param - 0.75) * 4;
 339                         }
 340                         break;
 341
 342                 case HAACK:
 343                         int1 = vonKarmanInterpolator;
 344                         int2 = lvHaackInterpolator;
 345                         p = param * 3;
 346                         break;
 347
 348                 default:
 349                         throw new UnsupportedOperationException("Unknown transition shape: " + shape);
 350                 }
 351
 352                 if (p < 0 || p > 1.00001) {
 353                         throw new BugException("Inconsistent parameter value p=" + p + " shape=" + shape);
 354                 }
 355
 356
 357                 // Check for parameterized shape and interpolate if necessary
 358                 if (int2 != null) {
 359                         LinearInterpolator int3 = new LinearInterpolator();
 360                         for (double m : int1.getXPoints()) {
 361                                 int3.addPoint(m, p * int2.getValue(m) + (1 - p) * int1.getValue(m));
 362                         }
 363                         for (double m : int2.getXPoints()) {
 364                                 int3.addPoint(m, p * int2.getValue(m) + (1 - p) * int1.getValue(m));
 365                         }
 366                         int1 = int3;
 367                 }
 368
 369                 // Extrapolate for fineness ratio if necessary
 370                 if (int1 != null) {
 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));
 375                         }
 376                 }
 377
 378
 379                 /*
 380                  * Now the transonic/supersonic region is ok.  We still need to interpolate
 381                  * the subsonic region, if the values are non-zero.
 382                  */
 383
 384                 double min = interpolator.getXPoints()[0];
 385                 double minValue = interpolator.getValue(min);
 386                 if (minValue < 0.001) {
 387                         // No interpolation necessary
 388                         return;
 389                 }
 390
 391                 double cdMach0 = 0.8 * pow2(sinphi);
 392                 double minDeriv = (interpolator.getValue(min + 0.01) - minValue) / 0.01;
 393
 394                 // These should not occur, but might cause havoc for the interpolation
 395                 if ((cdMach0 >= minValue - 0.01) || (minDeriv <= 0.01)) {
 396                         return;
 397                 }
 398
 399                 // Cd = a*M^b + cdMach0
 400                 double a = minValue - cdMach0;
 401                 double b = minDeriv / a;
 402
 403                 for (double m = 0; m < minValue; m += 0.05) {
 404                         interpolator.addPoint(m, a * Math.pow(m, b) + cdMach0);
 405                 }
 406         }
 407
 408
 409         private static final PolyInterpolator conicalPolyInterpolator =
 410                         new PolyInterpolator(new double[] { 1.0, 1.3 }, new double[] { 1.0, 1.3 });
 411
 412         private static LinearInterpolator calculateOgiveNoseInterpolator(double param,
 413                         double sinphi) {
 414                 LinearInterpolator interpolator = new LinearInterpolator();
 415
 416                 // In the range M = 1 ... 1.3 use polynomial approximation
 417                 double cdMach1 = 2.1 * pow2(sinphi) + 0.6019 * sinphi;
 418
 419                 double[] poly = conicalPolyInterpolator.interpolator(
 420                                 1.0 * sinphi, cdMach1,
 421                                 4 / (GAMMA + 1) * (1 - 0.5 * cdMach1), -1.1341 * sinphi
 422                                 );
 423
 424                 // Shape parameter multiplier
 425                 double mul = 0.72 * pow2(param - 0.5) + 0.82;
 426
 427                 for (double m = 1; m < 1.3001; m += 0.02) {
 428                         interpolator.addPoint(m, mul * PolyInterpolator.eval(m, poly));
 429                 }
 430
 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)));
 434                 }
 435
 436                 return interpolator;
 437         }
 438
 439
 440
 441 }