1 package net.sf.openrocket.aerodynamics.barrowman;
3 import static java.lang.Math.*;
4 import static net.sf.openrocket.util.MathUtil.pow2;
6 import java.util.Arrays;
8 import net.sf.openrocket.aerodynamics.AerodynamicForces;
9 import net.sf.openrocket.aerodynamics.FlightConditions;
10 import net.sf.openrocket.aerodynamics.Warning;
11 import net.sf.openrocket.aerodynamics.WarningSet;
12 import net.sf.openrocket.rocketcomponent.FinSet;
13 import net.sf.openrocket.rocketcomponent.RocketComponent;
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;
20 public class FinSetCalc extends RocketComponentCalc {
21 private static final double STALL_ANGLE = (20 * Math.PI/180);
23 /** Number of divisions in the fin chords. */
24 protected static final int DIVISIONS = 48;
29 private FinSet component;
32 protected double macLength = Double.NaN; // MAC length
33 protected double macLead = Double.NaN; // MAC leading edge position
34 protected double macSpan = Double.NaN; // MAC spanwise position
35 protected double finArea = Double.NaN; // Fin area
36 protected double ar = Double.NaN; // Fin aspect ratio
37 protected double span = Double.NaN; // Fin span
38 protected double cosGamma = Double.NaN; // Cosine of midchord sweep angle
39 protected double cosGammaLead = Double.NaN; // Cosine of leading edge sweep angle
40 protected double rollSum = Double.NaN; // Roll damping sum term
42 protected double[] chordLead = new double[DIVISIONS];
43 protected double[] chordTrail = new double[DIVISIONS];
44 protected double[] chordLength = new double[DIVISIONS];
46 protected final WarningSet geometryWarnings = new WarningSet();
48 private double[] poly = new double[6];
51 public FinSetCalc(RocketComponent component) {
53 if (!(component instanceof FinSet)) {
54 throw new IllegalArgumentException("Illegal component type "+component);
56 this.component = (FinSet) component;
61 * Calculates the non-axial forces produced by the fins (normal and side forces,
62 * pitch, yaw and roll moments, CP position, CNa).
65 public void calculateNonaxialForces(FlightConditions conditions,
66 AerodynamicForces forces, WarningSet warnings) {
68 // Compute and cache the fin geometry
69 if (Double.isNaN(macLength)) {
70 calculateFinGeometry();
78 forces.cp = Coordinate.NUL;
81 forces.CrollForce = 0;
88 // Add warnings (radius/2 == diameter/4)
89 if (component.getThickness() > component.getBodyRadius()/2) {
90 warnings.add(Warning.THICK_FIN);
92 warnings.addAll(geometryWarnings);
96 //////// Calculate CNa. /////////
98 // One fin without interference (both sub- and supersonic):
99 double cna1 = calculateFinCNa1(conditions);
103 // Multiple fins with fin-fin interference
106 // TODO: MEDIUM: Take into account multiple fin sets
107 int fins = component.getFinCount();
108 double theta = conditions.getTheta();
109 double angle = component.getBaseRotation();
116 for (int i=0; i < fins; i++) {
117 mul += MathUtil.pow2(Math.sin(theta - angle));
118 angle += 2 * Math.PI / fins;
124 // multiplier 1.5, sinusoidal reduction of 15%
125 cna = cna1 * 1.5 * (1 - 0.15*pow2(Math.cos(1.5 * (theta-angle))));
129 // multiplier 2.0, sinusoidal reduction of 6%
130 cna = cna1 * 2.0 * (1 - 0.06*pow2(Math.sin(2 * (theta-angle))));
150 // Assume N/2 * 3/4 efficiency for more fins
151 cna = cna1 * fins * 3.0/8.0;
156 // Body-fin interference effect
157 double r = component.getBodyRadius();
158 double tau = r / (span+r);
159 if (Double.isNaN(tau) || Double.isInfinite(tau))
161 cna *= 1 + tau; // Classical Barrowman
162 // cna *= pow2(1 + tau); // Barrowman thesis (too optimistic??)
166 // TODO: LOW: check for fin tip mach cone interference
167 // (Barrowman thesis pdf-page 40)
169 // TODO: LOW: fin-fin mach cone effect, MIL-HDBK page 5-25
173 // Calculate CP position
174 double x = macLead + calculateCPPos(conditions) * macLength;
178 // Calculate roll forces, reduce forcing above stall angle
180 // Without body-fin interference effect:
181 // forces.CrollForce = fins * (macSpan+r) * cna1 * component.getCantAngle() /
182 // conditions.getRefLength();
183 // With body-fin interference effect:
184 forces.CrollForce = fins * (macSpan+r) * cna1 * (1+tau) * component.getCantAngle() /
185 conditions.getRefLength();
190 if (conditions.getAOA() > STALL_ANGLE) {
191 // System.out.println("Fin stalling in roll");
192 forces.CrollForce *= MathUtil.clamp(
193 1-(conditions.getAOA() - STALL_ANGLE)/(STALL_ANGLE/2), 0, 1);
195 forces.CrollDamp = calculateDampingMoment(conditions);
196 forces.Croll = forces.CrollForce - forces.CrollDamp;
200 // System.out.printf(component.getName() + ": roll rate:%.3f force:%.3f damp:%.3f " +
202 // conditions.getRollRate(), forces.CrollForce, forces.CrollDamp, forces.Croll);
205 forces.CN = cna * MathUtil.min(conditions.getAOA(), STALL_ANGLE);
206 forces.cp = new Coordinate(x, 0, 0, cna);
207 forces.Cm = forces.CN * x / conditions.getRefLength();
210 forces.Cside = cna1 * Math.cos(theta-angle) * Math.sin(theta-angle);
211 forces.Cyaw = forces.Cside * x / conditions.getRefLength();
221 * Returns the MAC length of the fin. This is required in the friction drag
224 * @return the MAC length of the fin.
226 public double getMACLength() {
227 // Compute and cache the fin geometry
228 if (Double.isNaN(macLength)) {
229 calculateFinGeometry();
236 public double getMidchordPos() {
237 // Compute and cache the fin geometry
238 if (Double.isNaN(macLength)) {
239 calculateFinGeometry();
243 return macLead + 0.5 * macLength;
249 * Pre-calculates the fin geometry values.
251 protected void calculateFinGeometry() {
253 span = component.getSpan();
254 finArea = component.getFinArea();
255 ar = 2 * pow2(span) / finArea;
257 Coordinate[] points = component.getFinPoints();
259 // Check for jagged edges
260 geometryWarnings.clear();
261 boolean down = false;
262 for (int i=1; i < points.length; i++) {
263 if ((points[i].y > points[i-1].y + 0.001) && down) {
264 geometryWarnings.add(Warning.JAGGED_EDGED_FIN);
267 if (points[i].y < points[i-1].y - 0.001) {
273 // Calculate the chord lead and trail positions and length
275 Arrays.fill(chordLead, Double.POSITIVE_INFINITY);
276 Arrays.fill(chordTrail, Double.NEGATIVE_INFINITY);
277 Arrays.fill(chordLength, 0);
279 for (int point=1; point < points.length; point++) {
280 double x1 = points[point-1].x;
281 double y1 = points[point-1].y;
282 double x2 = points[point].x;
283 double y2 = points[point].y;
285 if (MathUtil.equals(y1, y2))
288 int i1 = (int)(y1*1.0001/span*(DIVISIONS-1));
289 int i2 = (int)(y2*1.0001/span*(DIVISIONS-1));
290 i1 = MathUtil.clamp(i1, 0, DIVISIONS-1);
291 i2 = MathUtil.clamp(i2, 0, DIVISIONS-1);
298 for (int i = i1; i <= i2; i++) {
299 // Intersection point (x,y)
300 double y = i*span/(DIVISIONS-1);
301 double x = (y-y2)/(y1-y2)*x1 + (y1-y)/(y1-y2)*x2;
302 if (x < chordLead[i])
304 if (x > chordTrail[i])
307 // TODO: LOW: If fin point exactly on chord line, might be counted twice:
316 // Check and correct any inconsistencies
317 for (int i=0; i < DIVISIONS; i++) {
318 if (Double.isInfinite(chordLead[i]) || Double.isInfinite(chordTrail[i]) ||
319 Double.isNaN(chordLead[i]) || Double.isNaN(chordTrail[i])) {
323 if (chordLength[i] < 0 || Double.isNaN(chordLength[i])) {
326 if (chordLength[i] > chordTrail[i] - chordLead[i]) {
327 chordLength[i] = chordTrail[i] - chordLead[i];
332 /* Calculate fin properties:
334 * macLength // MAC length
335 * macLead // MAC leading edge position
336 * macSpan // MAC spanwise position
337 * ar // Fin aspect ratio (already set)
338 * span // Fin span (already set)
347 double radius = component.getBodyRadius();
349 final double dy = span/(DIVISIONS-1);
350 for (int i=0; i < DIVISIONS; i++) {
351 double length = chordTrail[i] - chordLead[i];
354 macLength += length * length;
355 macSpan += y * length;
356 macLead += chordLead[i] * length;
358 rollSum += chordLength[i] * pow2(radius + y);
361 double dx = (chordTrail[i]+chordLead[i])/2 - (chordTrail[i-1]+chordLead[i-1])/2;
362 cosGamma += dy/MathUtil.hypot(dx, dy);
364 dx = chordLead[i] - chordLead[i-1];
365 cosGammaLead += dy/MathUtil.hypot(dx, dy);
378 cosGamma /= (DIVISIONS-1);
379 cosGammaLead /= (DIVISIONS-1);
383 /////////////// CNa1 calculation ////////////////
385 private static final double CNA_SUBSONIC = 0.9;
386 private static final double CNA_SUPERSONIC = 1.5;
387 private static final double CNA_SUPERSONIC_B = pow(pow2(CNA_SUPERSONIC)-1, 1.5);
388 private static final double GAMMA = 1.4;
389 private static final LinearInterpolator K1, K2, K3;
390 private static final PolyInterpolator cnaInterpolator = new PolyInterpolator(
391 new double[] { CNA_SUBSONIC, CNA_SUPERSONIC },
392 new double[] { CNA_SUBSONIC, CNA_SUPERSONIC },
393 new double[] { CNA_SUBSONIC }
395 /* Pre-calculate the values for K1, K2 and K3 */
398 int n = (int)((5.0-CNA_SUPERSONIC)*10);
399 double[] x = new double[n];
400 double[] k1 = new double[n];
401 double[] k2 = new double[n];
402 double[] k3 = new double[n];
403 for (int i=0; i<n; i++) {
404 double M = CNA_SUPERSONIC + i*0.1;
405 double beta = sqrt(M*M - 1);
408 k2[i] = ((GAMMA+1)*pow(M, 4) - 4*pow2(beta)) / (4*pow(beta,4));
409 k3[i] = ((GAMMA+1)*pow(M, 8) + (2*pow2(GAMMA) - 7*GAMMA - 5) * pow(M,6) +
410 10*(GAMMA+1)*pow(M,4) + 8) / (6*pow(beta,7));
412 K1 = new LinearInterpolator(x,k1);
413 K2 = new LinearInterpolator(x,k2);
414 K3 = new LinearInterpolator(x,k3);
416 // System.out.println("K1[m="+CNA_SUPERSONIC+"] = "+k1[0]);
417 // System.out.println("K2[m="+CNA_SUPERSONIC+"] = "+k2[0]);
418 // System.out.println("K3[m="+CNA_SUPERSONIC+"] = "+k3[0]);
422 protected double calculateFinCNa1(FlightConditions conditions) {
423 double mach = conditions.getMach();
424 double ref = conditions.getRefArea();
425 double alpha = MathUtil.min(conditions.getAOA(),
426 Math.PI - conditions.getAOA(), STALL_ANGLE);
429 if (mach <= CNA_SUBSONIC) {
430 return 2*Math.PI*pow2(span)/(1 + sqrt(1 + (1-pow2(mach))*
431 pow2(pow2(span)/(finArea*cosGamma)))) / ref;
435 if (mach >= CNA_SUPERSONIC) {
436 return finArea * (K1.getValue(mach) + K2.getValue(mach)*alpha +
437 K3.getValue(mach)*pow2(alpha)) / ref;
440 // Transonic case, interpolate
444 double sq = sqrt(1 + (1-pow2(CNA_SUBSONIC)) * pow2(span*span/(finArea*cosGamma)));
445 subV = 2*Math.PI*pow2(span)/ref / (1+sq);
446 subD = 2*mach*Math.PI*pow(span,6) / (pow2(finArea*cosGamma) * ref *
449 superV = finArea * (K1.getValue(CNA_SUPERSONIC) + K2.getValue(CNA_SUPERSONIC)*alpha +
450 K3.getValue(CNA_SUPERSONIC)*pow2(alpha)) / ref;
451 superD = -finArea/ref * 2*CNA_SUPERSONIC / CNA_SUPERSONIC_B;
453 // System.out.println("subV="+subV+" superV="+superV+" subD="+subD+" superD="+superD);
455 return cnaInterpolator.interpolate(mach, subV, superV, subD, superD, 0);
461 private double calculateDampingMoment(FlightConditions conditions) {
462 double rollRate = conditions.getRollRate();
464 if (Math.abs(rollRate) < 0.1)
467 double mach = conditions.getMach();
468 double radius = component.getBodyRadius();
469 double absRate = Math.abs(rollRate);
473 * At low speeds and relatively large roll rates (i.e. near apogee) the
474 * fin tips rotate well above stall angle. In this case sum the chords
477 if (absRate * (radius + span) / conditions.getVelocity() > 15*Math.PI/180) {
479 for (int i=0; i < DIVISIONS; i++) {
480 double dist = radius + span*i/DIVISIONS;
481 double aoa = Math.min(absRate*dist/conditions.getVelocity(), 15*Math.PI/180);
482 sum += chordLength[i] * dist * aoa;
484 sum = sum * (span/DIVISIONS) * 2*Math.PI/conditions.getBeta() /
485 (conditions.getRefArea() * conditions.getRefLength());
487 // System.out.println("SPECIAL: " +
488 // (MathUtil.sign(rollRate) *component.getFinCount() * sum));
489 return MathUtil.sign(rollRate) * component.getFinCount() * sum;
494 if (mach <= CNA_SUBSONIC) {
495 // System.out.println("BASIC: "+
496 // (component.getFinCount() * 2*Math.PI * rollRate * rollSum /
497 // (conditions.getRefArea() * conditions.getRefLength() *
498 // conditions.getVelocity() * conditions.getBeta())));
500 return component.getFinCount() * 2*Math.PI * rollRate * rollSum /
501 (conditions.getRefArea() * conditions.getRefLength() *
502 conditions.getVelocity() * conditions.getBeta());
504 if (mach >= CNA_SUPERSONIC) {
506 double vel = conditions.getVelocity();
507 double k1 = K1.getValue(mach);
508 double k2 = K2.getValue(mach);
509 double k3 = K3.getValue(mach);
513 for (int i=0; i < DIVISIONS; i++) {
514 double y = i*span/(DIVISIONS-1);
515 double angle = rollRate * (radius+y) / vel;
517 sum += (k1 * angle + k2 * angle*angle + k3 * angle*angle*angle)
518 * chordLength[i] * (radius+y);
521 return component.getFinCount() * sum * span/(DIVISIONS-1) /
522 (conditions.getRefArea() * conditions.getRefLength());
525 // Transonic, do linear interpolation
527 FlightConditions cond = conditions.clone();
528 cond.setMach(CNA_SUBSONIC - 0.01);
529 double subsonic = calculateDampingMoment(cond);
530 cond.setMach(CNA_SUPERSONIC + 0.01);
531 double supersonic = calculateDampingMoment(cond);
533 return subsonic * (CNA_SUPERSONIC - mach)/(CNA_SUPERSONIC - CNA_SUBSONIC) +
534 supersonic * (mach - CNA_SUBSONIC)/(CNA_SUPERSONIC - CNA_SUBSONIC);
541 * Return the relative position of the CP along the mean aerodynamic chord.
542 * Below mach 0.5 it is at the quarter chord, above mach 2 calculated using an
543 * empirical formula, between these two using an interpolation polynomial.
545 * @param cond Mach speed used
546 * @return CP position along the MAC
548 private double calculateCPPos(FlightConditions cond) {
549 double m = cond.getMach();
551 // At subsonic speeds CP at quarter chord
555 // At supersonic speeds use empirical formula
556 double beta = cond.getBeta();
557 return (ar * beta - 0.67) / (2*ar*beta - 1);
560 // In between use interpolation polynomial
564 for (int i=0; i < poly.length; i++) {
573 * Calculate CP position interpolation polynomial coefficients from the
574 * fin geometry. This is a fifth order polynomial that satisfies
583 * where f(M) = (ar*sqrt(M^2-1) - 0.67) / (2*ar*sqrt(M^2-1) - 1).
585 * The values were calculated analytically in Mathematica. The coefficients
586 * are used as poly[0] + poly[1]*x + poly[2]*x^2 + ...
588 private void calculatePoly() {
589 double denom = pow2(1 - 3.4641*ar); // common denominator
591 poly[5] = (-1.58025 * (-0.728769 + ar) * (-0.192105 + ar)) / denom;
592 poly[4] = (12.8395 * (-0.725688 + ar) * (-0.19292 + ar)) / denom;
593 poly[3] = (-39.5062 * (-0.72074 + ar) * (-0.194245 + ar)) / denom;
594 poly[2] = (55.3086 * (-0.711482 + ar) * (-0.196772 + ar)) / denom;
595 poly[1] = (-31.6049 * (-0.705375 + ar) * (-0.198476 + ar)) / denom;
596 poly[0] = (9.16049 * (-0.588838 + ar) * (-0.20624 + ar)) / denom;
600 // @SuppressWarnings("null")
601 // public static void main(String arg[]) {
602 // Rocket rocket = TestRocket.makeRocket();
603 // FinSet finset = null;
605 // Iterator<RocketComponent> iter = rocket.deepIterator();
606 // while (iter.hasNext()) {
607 // RocketComponent c = iter.next();
608 // if (c instanceof FinSet) {
609 // finset = (FinSet)c;
614 // ((TrapezoidFinSet)finset).setHeight(0.10);
615 // ((TrapezoidFinSet)finset).setRootChord(0.10);
616 // ((TrapezoidFinSet)finset).setTipChord(0.10);
617 // ((TrapezoidFinSet)finset).setSweep(0.0);
620 // FinSetCalc calc = new FinSetCalc(finset);
622 // calc.calculateFinGeometry();
623 // FlightConditions cond = new FlightConditions(new Configuration(rocket));
624 // for (double m=0; m < 3; m+=0.05) {
626 // cond.setAOA(0.0*Math.PI/180);
627 // double cna = calc.calculateFinCNa1(cond);
628 // System.out.printf("%5.2f "+cna+"\n", m);
635 public double calculatePressureDragForce(FlightConditions conditions,
636 double stagnationCD, double baseCD, WarningSet warnings) {
638 // Compute and cache the fin geometry
639 if (Double.isNaN(cosGammaLead)) {
640 calculateFinGeometry();
645 FinSet.CrossSection profile = component.getCrossSection();
646 double mach = conditions.getMach();
649 // Pressure fore-drag
650 if (profile == FinSet.CrossSection.AIRFOIL ||
651 profile == FinSet.CrossSection.ROUNDED) {
653 // Round leading edge
655 drag = Math.pow(1 - pow2(mach), -0.417) - 1;
656 } else if (mach < 1) {
657 drag = 1 - 1.785 * (mach-0.9);
659 drag = 1.214 - 0.502/pow2(mach) + 0.1095/pow2(pow2(mach));
662 } else if (profile == FinSet.CrossSection.SQUARE) {
665 throw new UnsupportedOperationException("Unsupported fin profile: "+profile);
668 // Slanted leading edge
669 drag *= pow2(cosGammaLead);
671 // Trailing edge drag
672 if (profile == FinSet.CrossSection.SQUARE) {
674 } else if (profile == FinSet.CrossSection.ROUNDED) {
677 // Airfoil assumed to have zero base drag
680 // Scale to correct reference area
681 drag *= component.getFinCount() * component.getSpan() * component.getThickness() /
682 conditions.getRefArea();