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