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

   1 package net.sf.openrocket.aerodynamics.barrowman;
   2
   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.BugException;
  14 import net.sf.openrocket.util.Coordinate;
  15 import net.sf.openrocket.util.LinearInterpolator;
  16 import net.sf.openrocket.util.MathUtil;
  17 import net.sf.openrocket.util.PolyInterpolator;
  18
  19
  20
  21 /**
  22  * Calculates the aerodynamic properties of a <code>SymmetricComponent</code>.
  23  * <p>
  24  * CP and CNa are calculated by the Barrowman method extended to account for body lift
  25  * by the method presented by Galejs.  Supersonic CNa and CP are assumed to be the
  26  * same as the subsonic values.
  27  *
  28  *
  29  * @author Sampo Niskanen <sampo.niskanen@iki.fi>
  30  */
  31 public class SymmetricComponentCalc extends RocketComponentCalc {
  32
  33         public static final double BODY_LIFT_K = 1.1;
  34
  35         private final SymmetricComponent component;
  36
  37         private final double length;
  38         private final double r1, r2;
  39         private final double fineness;
  40         private final Transition.Shape shape;
  41         private final double param;
  42         private final double area;
  43
  44         public SymmetricComponentCalc(RocketComponent c) {
  45                 super(c);
  46                 if (!(c instanceof SymmetricComponent)) {
  47                         throw new IllegalArgumentException("Illegal component type "+c);
  48                 }
  49                 this.component = (SymmetricComponent) c;
  50
  51
  52                 length = component.getLength();
  53                 r1 = component.getForeRadius();
  54                 r2 = component.getAftRadius();
  55
  56                 fineness = length / (2*Math.abs(r2-r1));
  57
  58                 if (component instanceof BodyTube) {
  59                         shape = null;
  60                         param = 0;
  61                         area = 0;
  62                 } else if (component instanceof Transition) {
  63                         shape = ((Transition)component).getType();
  64                         param = ((Transition)component).getShapeParameter();
  65                         area = Math.abs(Math.PI * (r1*r1 - r2*r2));
  66                 } else {
  67                         throw new UnsupportedOperationException("Unknown component type " +
  68                                         component.getComponentName());
  69                 }
  70         }
  71
  72
  73         private boolean isTube = false;
  74         private double cnaCache = Double.NaN;
  75         private double cpCache = Double.NaN;
  76
  77
  78         /**
  79          * Calculates the non-axial forces produced by the fins (normal and side forces,
  80          * pitch, yaw and roll moments, CP position, CNa).
  81          * <p>
  82          * This method uses the Barrowman method for CP and CNa calculation and the
  83          * extension presented by Galejs for the effect of body lift.
  84          * <p>
  85          * The CP and CNa at supersonic speeds are assumed to be the same as those at
  86          * subsonic speeds.
  87          */
  88         @Override
  89         public void calculateNonaxialForces(FlightConditions conditions,
  90                         AerodynamicForces forces, WarningSet warnings) {
  91
  92                 // Pre-calculate and store the results
  93                 if (Double.isNaN(cnaCache)) {
  94                         final double r0 = component.getForeRadius();
  95                         final double r1 = component.getAftRadius();
  96
  97                         if (MathUtil.equals(r0, r1)) {
  98                                 isTube = true;
  99                                 cnaCache = 0;
 100                         } else {
 101                                 isTube = false;
 102
 103                                 final double A0 = Math.PI * pow2(r0);
 104                                 final double A1 = Math.PI * pow2(r1);
 105
 106                                 cnaCache = 2 * (A1 - A0);
 107                                 System.out.println("cnaCache = "+cnaCache);
 108                                 cpCache = (component.getLength() * A1 - component.getFullVolume()) / (A1 - A0);
 109                         }
 110                 }
 111
 112                 Coordinate cp;
 113
 114                 // If fore == aft, only body lift is encountered
 115                 if (isTube) {
 116                         cp = getLiftCP(conditions, warnings);
 117                 } else {
 118                         cp = new Coordinate(cpCache,0,0,cnaCache * conditions.getSincAOA() /
 119                                         conditions.getRefArea()).average(getLiftCP(conditions,warnings));
 120                 }
 121
 122                 forces.cp = cp;
 123                 forces.CNa = cp.weight;
 124                 forces.CN = forces.CNa * conditions.getAOA();
 125                 forces.Cm = forces.CN * cp.x / conditions.getRefLength();
 126                 forces.Croll = 0;
 127                 forces.CrollDamp = 0;
 128                 forces.CrollForce = 0;
 129                 forces.Cside = 0;
 130                 forces.Cyaw = 0;
 131
 132
 133                 // Add warning on supersonic flight
 134                 if (conditions.getMach() > 1.1) {
 135                         warnings.add("Body calculations may not be entirely accurate at supersonic speeds.");
 136                 }
 137
 138         }
 139
 140
 141
 142         /**
 143          * Calculate the body lift effect according to Galejs.
 144          */
 145         protected Coordinate getLiftCP(FlightConditions conditions, WarningSet warnings) {
 146                 double area = component.getComponentPlanformArea();
 147                 double center = component.getComponentPlanformCenter();
 148
 149                 /*
 150                  * Without this extra multiplier the rocket may become unstable at apogee
 151                  * when turning around, and begin oscillating horizontally.  During the flight
 152                  * of the rocket this has no effect.  It is effective only when AOA > 45 deg
 153                  * and the velocity is less than 15 m/s.
 154                  */
 155                 double mul = 1;
 156                 if ((conditions.getMach() < 0.05) && (conditions.getAOA() > Math.PI/4)) {
 157                         mul = pow2(conditions.getMach() / 0.05);
 158                 }
 159
 160                 return new Coordinate(center, 0, 0, mul*BODY_LIFT_K * area/conditions.getRefArea() *
 161                                 conditions.getSinAOA() * conditions.getSincAOA());  // sin(aoa)^2 / aoa
 162         }
 163
 164
 165
 166         private LinearInterpolator interpolator = null;
 167
 168         @Override
 169         public double calculatePressureDragForce(FlightConditions conditions,
 170                         double stagnationCD, double baseCD, WarningSet warnings) {
 171
 172                 if (component instanceof BodyTube)
 173                         return 0;
 174
 175                 if (!(component instanceof Transition)) {
 176                         throw new BugException("Pressure calculation of unknown type: "+
 177                                         component.getComponentName());
 178                 }
 179
 180                 // Check for simple cases first
 181                 if (r1 == r2)
 182                         return 0;
 183
 184                 if (length < 0.001) {
 185                         if (r1 < r2) {
 186                                 return stagnationCD * area / conditions.getRefArea();
 187                         } else {
 188                                 return baseCD * area / conditions.getRefArea();
 189                         }
 190                 }
 191
 192
 193                 // Boattail drag computed directly from base drag
 194                 if (r2 < r1) {
 195                         if (fineness >= 3)
 196                                 return 0;
 197                         double cd = baseCD * area / conditions.getRefArea();
 198                         if (fineness <= 1)
 199                                 return cd;
 200                         return cd * (3-fineness)/2;
 201                 }
 202
 203
 204                 assert(r1 < r2);  // Tube and boattail have been checked already
 205
 206
 207                 // All nose cones and shoulders from pre-calculated and interpolating
 208                 if (interpolator == null) {
 209                         calculateNoseInterpolator();
 210                 }
 211
 212                 return interpolator.getValue(conditions.getMach()) * area / conditions.getRefArea();
 213         }
 214
 215
 216
 217         /*
 218          * Experimental values of pressure drag for different nose cone shapes with a fineness
 219          * ratio of 3.  The data is taken from 'Collection of Zero-Lift Drag Data on Bodies
 220          * of Revolution from Free-Flight Investigations', NASA TR-R-100, NTRS 19630004995,
 221          * page 16.
 222          *
 223          * This data is extrapolated for other fineness ratios.
 224          */
 225
 226         private static final LinearInterpolator ellipsoidInterpolator = new LinearInterpolator(
 227                         new double[] {  1.2,  1.25,   1.3,   1.4,   1.6,   2.0,   2.4 },
 228                         new double[] {0.110, 0.128, 0.140, 0.148, 0.152, 0.159, 0.162 /* constant */ }
 229         );
 230         private static final LinearInterpolator x14Interpolator = new LinearInterpolator(
 231                         new double[] {  1.2,   1.3,   1.4,   1.6,   1.8,   2.2,   2.6,   3.0,   3.6},
 232                         new double[] {0.140, 0.156, 0.169, 0.192, 0.206, 0.227, 0.241, 0.249, 0.252}
 233         );
 234         private static final LinearInterpolator x12Interpolator = new LinearInterpolator(
 235                         new double[] {0.925,  0.95,   1.0,  1.05,   1.1,   1.2,   1.3,   1.7,   2.0},
 236                         new double[] {    0, 0.014, 0.050, 0.060, 0.059, 0.081, 0.084, 0.085, 0.078}
 237         );
 238         private static final LinearInterpolator x34Interpolator = new LinearInterpolator(
 239                         new double[] { 0.8,   0.9,   1.0,  1.06,   1.2,   1.4,   1.6,   2.0,   2.8,   3.4},
 240                         new double[] {   0, 0.015, 0.078, 0.121, 0.110, 0.098, 0.090, 0.084, 0.078, 0.074}
 241         );
 242         private static final LinearInterpolator vonKarmanInterpolator = new LinearInterpolator(
 243                         new double[] { 0.9,  0.95,   1.0,  1.05,   1.1,   1.2,   1.4,   1.6,   2.0,   3.0},
 244                         new double[] {   0, 0.010, 0.027, 0.055, 0.070, 0.081, 0.095, 0.097, 0.091, 0.083}
 245         );
 246         private static final LinearInterpolator lvHaackInterpolator = new LinearInterpolator(
 247                         new double[] { 0.9,  0.95,   1.0,  1.05,   1.1,   1.2,   1.4,   1.6,   2.0 },
 248                         new double[] {   0, 0.010, 0.024, 0.066, 0.084, 0.100, 0.114, 0.117, 0.113 }
 249         );
 250         private static final LinearInterpolator parabolicInterpolator = new LinearInterpolator(
 251                         new double[] {0.95, 0.975,   1.0,  1.05,   1.1,   1.2,   1.4,   1.7},
 252                         new double[] {   0, 0.016, 0.041, 0.092, 0.109, 0.119, 0.113, 0.108}
 253         );
 254         private static final LinearInterpolator parabolic12Interpolator = new LinearInterpolator(
 255                         new double[] { 0.8,   0.9,  0.95,   1.0,  1.05,   1.1,   1.3,   1.5,   1.8},
 256                         new double[] {   0, 0.016, 0.042, 0.100, 0.126, 0.125, 0.100, 0.090, 0.088}
 257         );
 258         private static final LinearInterpolator parabolic34Interpolator = new LinearInterpolator(
 259                         new double[] { 0.9,  0.95,   1.0,  1.05,   1.1,   1.2,   1.4,   1.7},
 260                         new double[] {   0, 0.023, 0.073, 0.098, 0.107, 0.106, 0.089, 0.082}
 261         );
 262         private static final LinearInterpolator bluntInterpolator = new LinearInterpolator();
 263         static {
 264                 for (double m=0; m<3; m+=0.05)
 265                         bluntInterpolator.addPoint(m, BarrowmanCalculator.calculateStagnationCD(m));
 266         }
 267
 268         /**
 269          * Calculate the LinearInterpolator 'interpolator'.  After this call, if can be used
 270          * to get the pressure drag coefficient at any Mach number.
 271          *
 272          * First, the transonic/supersonic region is computed.  For conical and ogive shapes
 273          * this is calculated directly.  For other shapes, the values for fineness-ratio 3
 274          * transitions are taken from the experimental values stored above (for parameterized
 275          * shapes the values are interpolated between the parameter values).  These are then
 276          * extrapolated to the current fineness ratio.
 277          *
 278          * Finally, if the first data points in the interpolator are not zero, the subsonic
 279          * region is interpolated in the form   Cd = a*M^b + Cd(M=0).
 280          */
 281         @SuppressWarnings("null")
 282         private void calculateNoseInterpolator() {
 283                 LinearInterpolator int1=null, int2=null;
 284                 double p = 0;
 285
 286                 interpolator = new LinearInterpolator();
 287
 288                 double r = component.getRadius(0.99*length);
 289                 double sinphi = (r2-r)/MathUtil.hypot(r2-r, 0.01*length);
 290
 291                 /*
 292                  * Take into account nose cone shape.  Conical and ogive generate the interpolator
 293                  * directly.  Others store a interpolator for fineness ratio 3 into int1, or
 294                  * for parameterized shapes store the bounding fineness ratio 3 interpolators into
 295                  * int1 and int2 and set 0 <= p <= 1 according to the bounds.
 296                  */
 297                 switch (shape) {
 298                 case CONICAL:
 299                         interpolator = calculateOgiveNoseInterpolator(0, sinphi);  // param==0 -> conical
 300                         break;
 301
 302                 case OGIVE:
 303                         interpolator = calculateOgiveNoseInterpolator(param, sinphi);
 304                         break;
 305
 306                 case ELLIPSOID:
 307                         int1 = ellipsoidInterpolator;
 308                         break;
 309
 310                 case POWER:
 311                         if (param <= 0.25) {
 312                                 int1 = bluntInterpolator;
 313                                 int2 = x14Interpolator;
 314                                 p = param*4;
 315                         } else if (param <= 0.5) {
 316                                 int1 = x14Interpolator;
 317                                 int2 = x12Interpolator;
 318                                 p = (param-0.25)*4;
 319                         } else if (param <= 0.75) {
 320                                 int1 = x12Interpolator;
 321                                 int2 = x34Interpolator;
 322                                 p = (param-0.5)*4;
 323                         } else {
 324                                 int1 = x34Interpolator;
 325                                 int2 = calculateOgiveNoseInterpolator(0, 1/Math.sqrt(1+4*pow2(fineness)));
 326                                 p = (param-0.75)*4;
 327                         }
 328                         break;
 329
 330                 case PARABOLIC:
 331                         if (param <= 0.5) {
 332                                 int1 = calculateOgiveNoseInterpolator(0, 1/Math.sqrt(1+4*pow2(fineness)));
 333                                 int2 = parabolic12Interpolator;
 334                                 p = param*2;
 335                         } else if (param <= 0.75) {
 336                                 int1 = parabolic12Interpolator;
 337                                 int2 = parabolic34Interpolator;
 338                                 p = (param-0.5)*4;
 339                         } else {
 340                                 int1 = parabolic34Interpolator;
 341                                 int2 = parabolicInterpolator;
 342                                 p = (param-0.75)*4;
 343                         }
 344                         break;
 345
 346                 case HAACK:
 347                         int1 = vonKarmanInterpolator;
 348                         int2 = lvHaackInterpolator;
 349                         p = param*3;
 350                         break;
 351
 352                 default:
 353                         throw new UnsupportedOperationException("Unknown transition shape: "+shape);
 354                 }
 355
 356                 assert(p >= 0);
 357                 assert(p <= 1.001);
 358
 359
 360                 // Check for parameterized shape and interpolate if necessary
 361                 if (int2 != null) {
 362                         LinearInterpolator int3 = new LinearInterpolator();
 363                         for (double m: int1.getXPoints()) {
 364                                 int3.addPoint(m, p*int2.getValue(m) + (1-p)*int1.getValue(m));
 365                         }
 366                         for (double m: int2.getXPoints()) {
 367                                 int3.addPoint(m, p*int2.getValue(m) + (1-p)*int1.getValue(m));
 368                         }
 369                         int1 = int3;
 370                 }
 371
 372                 // Extrapolate for fineness ratio if necessary
 373                 if (int1 != null) {
 374                         double log4 = Math.log(fineness+1) / Math.log(4);
 375                         for (double m: int1.getXPoints()) {
 376                                 double stag = bluntInterpolator.getValue(m);
 377                                 interpolator.addPoint(m, stag*Math.pow(int1.getValue(m)/stag, log4));
 378                         }
 379                 }
 380
 381
 382                 /*
 383                  * Now the transonic/supersonic region is ok.  We still need to interpolate
 384                  * the subsonic region, if the values are non-zero.
 385                  */
 386
 387                 double min = interpolator.getXPoints()[0];
 388                 double minValue = interpolator.getValue(min);
 389                 if (minValue < 0.001) {
 390                         // No interpolation necessary
 391                         return;
 392                 }
 393
 394                 double cdMach0 = 0.8 * pow2(sinphi);
 395                 double minDeriv = (interpolator.getValue(min+0.01) - minValue)/0.01;
 396
 397                 // These should not occur, but might cause havoc for the interpolation
 398                 if ((cdMach0 >= minValue-0.01) || (minDeriv <= 0.01)) {
 399                         return;
 400                 }
 401
 402                 // Cd = a*M^b + cdMach0
 403                 double a = minValue - cdMach0;
 404                 double b = minDeriv / a;
 405
 406                 for (double m=0; m < minValue; m+= 0.05) {
 407                         interpolator.addPoint(m, a*Math.pow(m, b) + cdMach0);
 408                 }
 409         }
 410
 411
 412         private static final PolyInterpolator conicalPolyInterpolator =
 413                 new PolyInterpolator(new double[] {1.0, 1.3}, new double[] {1.0, 1.3});
 414
 415         private static LinearInterpolator calculateOgiveNoseInterpolator(double param,
 416                         double sinphi) {
 417                 LinearInterpolator interpolator = new LinearInterpolator();
 418
 419                 // In the range M = 1 ... 1.3 use polynomial approximation
 420                 double cdMach1 = 2.1*pow2(sinphi) + 0.6019*sinphi;
 421
 422                 double[] poly = conicalPolyInterpolator.interpolator(
 423                                 1.0*sinphi, cdMach1,
 424                                 4/(GAMMA+1) * (1 - 0.5*cdMach1), -1.1341*sinphi
 425                 );
 426
 427                 // Shape parameter multiplier
 428                 double mul = 0.72 * pow2(param-0.5) + 0.82;
 429
 430                 for (double m = 1; m < 1.3001; m += 0.02) {
 431                         interpolator.addPoint(m, mul * PolyInterpolator.eval(m, poly));
 432                 }
 433
 434                 // Above M = 1.3 use direct formula
 435                 for (double m = 1.32; m < 4; m += 0.02) {
 436                         interpolator.addPoint(m, mul * (2.1*pow2(sinphi) + 0.5*sinphi/Math.sqrt(m*m - 1)));
 437                 }
 438
 439                 return interpolator;
 440         }
 441
 442
 443
 444 }