1 package net.sf.openrocket.aerodynamics.barrowman;
3 import static java.lang.Math.pow;
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.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;
21 public class FinSetCalc extends RocketComponentCalc {
22 private static final double STALL_ANGLE = (20 * Math.PI / 180);
24 /** Number of divisions in the fin chords. */
25 protected static final int DIVISIONS = 48;
29 protected double macLength = Double.NaN; // MAC length
30 protected double macLead = Double.NaN; // MAC leading edge position
31 protected double macSpan = Double.NaN; // MAC spanwise position
32 protected double finArea = Double.NaN; // Fin area
33 protected double ar = Double.NaN; // Fin aspect ratio
34 protected double span = Double.NaN; // Fin span
35 protected double cosGamma = Double.NaN; // Cosine of midchord sweep angle
36 protected double cosGammaLead = Double.NaN; // Cosine of leading edge sweep angle
37 protected double rollSum = Double.NaN; // Roll damping sum term
39 protected int interferenceFinCount = -1; // No. of fins in interference
41 protected double[] chordLead = new double[DIVISIONS];
42 protected double[] chordTrail = new double[DIVISIONS];
43 protected double[] chordLength = new double[DIVISIONS];
46 protected final WarningSet geometryWarnings = new WarningSet();
48 private double[] poly = new double[6];
50 private final double thickness;
51 private final double bodyRadius;
52 private final int finCount;
53 private final double baseRotation;
54 private final double cantAngle;
55 private final FinSet.CrossSection crossSection;
57 public FinSetCalc(RocketComponent component) {
59 if (!(component instanceof FinSet)) {
60 throw new IllegalArgumentException("Illegal component type " + component);
63 FinSet fin = (FinSet) component;
64 thickness = fin.getThickness();
65 bodyRadius = fin.getBodyRadius();
66 finCount = fin.getFinCount();
67 baseRotation = fin.getBaseRotation();
68 cantAngle = fin.getCantAngle();
70 finArea = fin.getFinArea();
71 crossSection = fin.getCrossSection();
73 calculateFinGeometry(fin);
75 calculateInterferenceFinCount(fin);
80 * Calculates the non-axial forces produced by the fins (normal and side forces,
81 * pitch, yaw and roll moments, CP position, CNa).
84 public void calculateNonaxialForces(FlightConditions conditions,
85 AerodynamicForces forces, WarningSet warnings) {
92 forces.setCP(Coordinate.NUL);
94 forces.setCrollDamp(0);
95 forces.setCrollForce(0);
102 // Add warnings (radius/2 == diameter/4)
103 if (thickness > bodyRadius / 2) {
104 warnings.add(Warning.THICK_FIN);
106 warnings.addAll(geometryWarnings);
110 //////// Calculate CNa. /////////
112 // One fin without interference (both sub- and supersonic):
113 double cna1 = calculateFinCNa1(conditions);
117 // Multiple fins with fin-fin interference
119 double theta = conditions.getTheta();
120 double angle = baseRotation;
123 // Compute basic CNa without interference effects
124 if (finCount == 1 || finCount == 2) {
125 // Basic CNa from geometry
127 for (int i = 0; i < finCount; i++) {
128 mul += MathUtil.pow2(Math.sin(theta - angle));
129 angle += 2 * Math.PI / finCount;
133 // Basic CNa assuming full efficiency
134 cna = cna1 * finCount / 2.0;
138 // Take into account fin-fin interference effects
139 switch (interferenceFinCount) {
144 // No interference effect
164 // Assume 75% efficiency
166 warnings.add(Warning.PARALLEL_FINS);
171 * Used in 0.9.5 and earlier. Takes into account rotation angle for three
172 * and four fins, does not take into account interference from other fin sets.
179 for (int i=0; i < fins; i++) {
180 mul += MathUtil.pow2(Math.sin(theta - angle));
181 angle += 2 * Math.PI / fins;
187 // multiplier 1.5, sinusoidal reduction of 15%
188 cna = cna1 * 1.5 * (1 - 0.15*pow2(Math.cos(1.5 * (theta-angle))));
192 // multiplier 2.0, sinusoidal reduction of 6%
193 cna = cna1 * 2.0 * (1 - 0.06*pow2(Math.sin(2 * (theta-angle))));
213 // Assume N/2 * 3/4 efficiency for more fins
214 cna = cna1 * fins * 3.0/8.0;
219 // Body-fin interference effect
220 double r = bodyRadius;
221 double tau = r / (span + r);
222 if (Double.isNaN(tau) || Double.isInfinite(tau))
224 cna *= 1 + tau; // Classical Barrowman
225 // cna *= pow2(1 + tau); // Barrowman thesis (too optimistic??)
229 // TODO: LOW: check for fin tip mach cone interference
230 // (Barrowman thesis pdf-page 40)
232 // TODO: LOW: fin-fin mach cone effect, MIL-HDBK page 5-25
236 // Calculate CP position
237 double x = macLead + calculateCPPos(conditions) * macLength;
241 // Calculate roll forces, reduce forcing above stall angle
243 // Without body-fin interference effect:
244 // forces.CrollForce = fins * (macSpan+r) * cna1 * component.getCantAngle() /
245 // conditions.getRefLength();
246 // With body-fin interference effect:
247 forces.setCrollForce(finCount * (macSpan + r) * cna1 * (1 + tau) * cantAngle / conditions.getRefLength());
252 if (conditions.getAOA() > STALL_ANGLE) {
253 // System.out.println("Fin stalling in roll");
254 forces.setCrollForce(forces.getCrollForce() * MathUtil.clamp(
255 1 - (conditions.getAOA() - STALL_ANGLE) / (STALL_ANGLE / 2), 0, 1));
257 forces.setCrollDamp(calculateDampingMoment(conditions));
258 forces.setCroll(forces.getCrollForce() - forces.getCrollDamp());
262 // System.out.printf(component.getName() + ": roll rate:%.3f force:%.3f damp:%.3f " +
264 // conditions.getRollRate(), forces.CrollForce, forces.CrollDamp, forces.Croll);
267 forces.setCN(cna * MathUtil.min(conditions.getAOA(), STALL_ANGLE));
268 forces.setCP(new Coordinate(x, 0, 0, cna));
269 forces.setCm(forces.getCN() * x / conditions.getRefLength());
272 * TODO: HIGH: Compute actual side force and yaw moment.
273 * This is not currently performed because it produces strange results for
274 * stable rockets that have two fins in the front part of the fuselage,
275 * where the rocket flies at an ever-increasing angle of attack. This may
276 * be due to incorrect computation of pitch/yaw damping moments.
278 // if (fins == 1 || fins == 2) {
279 // forces.Cside = fins * cna1 * Math.cos(theta-angle) * Math.sin(theta-angle);
280 // forces.Cyaw = fins * forces.Cside * x / conditions.getRefLength();
293 * Returns the MAC length of the fin. This is required in the friction drag
296 * @return the MAC length of the fin.
298 public double getMACLength() {
302 public double getMidchordPos() {
303 return macLead + 0.5 * macLength;
309 * Pre-calculates the fin geometry values.
311 protected void calculateFinGeometry(FinSet component) {
313 span = component.getSpan();
314 finArea = component.getFinArea();
315 ar = 2 * pow2(span) / finArea;
317 Coordinate[] points = component.getFinPoints();
319 // Check for jagged edges
320 geometryWarnings.clear();
321 boolean down = false;
322 for (int i = 1; i < points.length; i++) {
323 if ((points[i].y > points[i - 1].y + 0.001) && down) {
324 geometryWarnings.add(Warning.JAGGED_EDGED_FIN);
327 if (points[i].y < points[i - 1].y - 0.001) {
333 // Calculate the chord lead and trail positions and length
335 Arrays.fill(chordLead, Double.POSITIVE_INFINITY);
336 Arrays.fill(chordTrail, Double.NEGATIVE_INFINITY);
337 Arrays.fill(chordLength, 0);
339 for (int point = 1; point < points.length; point++) {
340 double x1 = points[point - 1].x;
341 double y1 = points[point - 1].y;
342 double x2 = points[point].x;
343 double y2 = points[point].y;
345 if (MathUtil.equals(y1, y2))
348 int i1 = (int) (y1 * 1.0001 / span * (DIVISIONS - 1));
349 int i2 = (int) (y2 * 1.0001 / span * (DIVISIONS - 1));
350 i1 = MathUtil.clamp(i1, 0, DIVISIONS - 1);
351 i2 = MathUtil.clamp(i2, 0, DIVISIONS - 1);
358 for (int i = i1; i <= i2; i++) {
359 // Intersection point (x,y)
360 double y = i * span / (DIVISIONS - 1);
361 double x = (y - y2) / (y1 - y2) * x1 + (y1 - y) / (y1 - y2) * x2;
362 if (x < chordLead[i])
364 if (x > chordTrail[i])
367 // TODO: LOW: If fin point exactly on chord line, might be counted twice:
376 // Check and correct any inconsistencies
377 for (int i = 0; i < DIVISIONS; i++) {
378 if (Double.isInfinite(chordLead[i]) || Double.isInfinite(chordTrail[i]) ||
379 Double.isNaN(chordLead[i]) || Double.isNaN(chordTrail[i])) {
383 if (chordLength[i] < 0 || Double.isNaN(chordLength[i])) {
386 if (chordLength[i] > chordTrail[i] - chordLead[i]) {
387 chordLength[i] = chordTrail[i] - chordLead[i];
392 /* Calculate fin properties:
394 * macLength // MAC length
395 * macLead // MAC leading edge position
396 * macSpan // MAC spanwise position
397 * ar // Fin aspect ratio (already set)
398 * span // Fin span (already set)
407 double radius = component.getBodyRadius();
409 final double dy = span / (DIVISIONS - 1);
410 for (int i = 0; i < DIVISIONS; i++) {
411 double length = chordTrail[i] - chordLead[i];
414 macLength += length * length;
415 macSpan += y * length;
416 macLead += chordLead[i] * length;
418 rollSum += chordLength[i] * pow2(radius + y);
421 double dx = (chordTrail[i] + chordLead[i]) / 2 - (chordTrail[i - 1] + chordLead[i - 1]) / 2;
422 cosGamma += dy / MathUtil.hypot(dx, dy);
424 dx = chordLead[i] - chordLead[i - 1];
425 cosGammaLead += dy / MathUtil.hypot(dx, dy);
438 cosGamma /= (DIVISIONS - 1);
439 cosGammaLead /= (DIVISIONS - 1);
443 /////////////// CNa1 calculation ////////////////
445 private static final double CNA_SUBSONIC = 0.9;
446 private static final double CNA_SUPERSONIC = 1.5;
447 private static final double CNA_SUPERSONIC_B = pow(pow2(CNA_SUPERSONIC) - 1, 1.5);
448 private static final double GAMMA = 1.4;
449 private static final LinearInterpolator K1, K2, K3;
450 private static final PolyInterpolator cnaInterpolator = new PolyInterpolator(
451 new double[] { CNA_SUBSONIC, CNA_SUPERSONIC },
452 new double[] { CNA_SUBSONIC, CNA_SUPERSONIC },
453 new double[] { CNA_SUBSONIC }
455 /* Pre-calculate the values for K1, K2 and K3 */
458 int n = (int) ((5.0 - CNA_SUPERSONIC) * 10);
459 double[] x = new double[n];
460 double[] k1 = new double[n];
461 double[] k2 = new double[n];
462 double[] k3 = new double[n];
463 for (int i = 0; i < n; i++) {
464 double M = CNA_SUPERSONIC + i * 0.1;
465 double beta = MathUtil.safeSqrt(M * M - 1);
468 k2[i] = ((GAMMA + 1) * pow(M, 4) - 4 * pow2(beta)) / (4 * pow(beta, 4));
469 k3[i] = ((GAMMA + 1) * pow(M, 8) + (2 * pow2(GAMMA) - 7 * GAMMA - 5) * pow(M, 6) +
470 10 * (GAMMA + 1) * pow(M, 4) + 8) / (6 * pow(beta, 7));
472 K1 = new LinearInterpolator(x, k1);
473 K2 = new LinearInterpolator(x, k2);
474 K3 = new LinearInterpolator(x, k3);
476 // System.out.println("K1[m="+CNA_SUPERSONIC+"] = "+k1[0]);
477 // System.out.println("K2[m="+CNA_SUPERSONIC+"] = "+k2[0]);
478 // System.out.println("K3[m="+CNA_SUPERSONIC+"] = "+k3[0]);
482 protected double calculateFinCNa1(FlightConditions conditions) {
483 double mach = conditions.getMach();
484 double ref = conditions.getRefArea();
485 double alpha = MathUtil.min(conditions.getAOA(),
486 Math.PI - conditions.getAOA(), STALL_ANGLE);
489 if (mach <= CNA_SUBSONIC) {
490 return 2 * Math.PI * pow2(span) / (1 + MathUtil.safeSqrt(1 + (1 - pow2(mach)) *
491 pow2(pow2(span) / (finArea * cosGamma)))) / ref;
495 if (mach >= CNA_SUPERSONIC) {
496 return finArea * (K1.getValue(mach) + K2.getValue(mach) * alpha +
497 K3.getValue(mach) * pow2(alpha)) / ref;
500 // Transonic case, interpolate
504 double sq = MathUtil.safeSqrt(1 + (1 - pow2(CNA_SUBSONIC)) * pow2(span * span / (finArea * cosGamma)));
505 subV = 2 * Math.PI * pow2(span) / ref / (1 + sq);
506 subD = 2 * mach * Math.PI * pow(span, 6) / (pow2(finArea * cosGamma) * ref *
509 superV = finArea * (K1.getValue(CNA_SUPERSONIC) + K2.getValue(CNA_SUPERSONIC) * alpha +
510 K3.getValue(CNA_SUPERSONIC) * pow2(alpha)) / ref;
511 superD = -finArea / ref * 2 * CNA_SUPERSONIC / CNA_SUPERSONIC_B;
513 // System.out.println("subV="+subV+" superV="+superV+" subD="+subD+" superD="+superD);
515 return cnaInterpolator.interpolate(mach, subV, superV, subD, superD, 0);
521 private double calculateDampingMoment(FlightConditions conditions) {
522 double rollRate = conditions.getRollRate();
524 if (Math.abs(rollRate) < 0.1)
527 double mach = conditions.getMach();
528 double absRate = Math.abs(rollRate);
532 * At low speeds and relatively large roll rates (i.e. near apogee) the
533 * fin tips rotate well above stall angle. In this case sum the chords
536 if (absRate * (bodyRadius + span) / conditions.getVelocity() > 15 * Math.PI / 180) {
538 for (int i = 0; i < DIVISIONS; i++) {
539 double dist = bodyRadius + span * i / DIVISIONS;
540 double aoa = Math.min(absRate * dist / conditions.getVelocity(), 15 * Math.PI / 180);
541 sum += chordLength[i] * dist * aoa;
543 sum = sum * (span / DIVISIONS) * 2 * Math.PI / conditions.getBeta() /
544 (conditions.getRefArea() * conditions.getRefLength());
546 // System.out.println("SPECIAL: " +
547 // (MathUtil.sign(rollRate) *component.getFinCount() * sum));
548 return MathUtil.sign(rollRate) * finCount * sum;
553 if (mach <= CNA_SUBSONIC) {
554 // System.out.println("BASIC: "+
555 // (component.getFinCount() * 2*Math.PI * rollRate * rollSum /
556 // (conditions.getRefArea() * conditions.getRefLength() *
557 // conditions.getVelocity() * conditions.getBeta())));
559 return finCount * 2 * Math.PI * rollRate * rollSum /
560 (conditions.getRefArea() * conditions.getRefLength() *
561 conditions.getVelocity() * conditions.getBeta());
563 if (mach >= CNA_SUPERSONIC) {
565 double vel = conditions.getVelocity();
566 double k1 = K1.getValue(mach);
567 double k2 = K2.getValue(mach);
568 double k3 = K3.getValue(mach);
572 for (int i = 0; i < DIVISIONS; i++) {
573 double y = i * span / (DIVISIONS - 1);
574 double angle = rollRate * (bodyRadius + y) / vel;
576 sum += (k1 * angle + k2 * angle * angle + k3 * angle * angle * angle)
577 * chordLength[i] * (bodyRadius + y);
580 return finCount * sum * span / (DIVISIONS - 1) /
581 (conditions.getRefArea() * conditions.getRefLength());
584 // Transonic, do linear interpolation
586 FlightConditions cond = conditions.clone();
587 cond.setMach(CNA_SUBSONIC - 0.01);
588 double subsonic = calculateDampingMoment(cond);
589 cond.setMach(CNA_SUPERSONIC + 0.01);
590 double supersonic = calculateDampingMoment(cond);
592 return subsonic * (CNA_SUPERSONIC - mach) / (CNA_SUPERSONIC - CNA_SUBSONIC) +
593 supersonic * (mach - CNA_SUBSONIC) / (CNA_SUPERSONIC - CNA_SUBSONIC);
600 * Return the relative position of the CP along the mean aerodynamic chord.
601 * Below mach 0.5 it is at the quarter chord, above mach 2 calculated using an
602 * empirical formula, between these two using an interpolation polynomial.
604 * @param cond Mach speed used
605 * @return CP position along the MAC
607 private double calculateCPPos(FlightConditions cond) {
608 double m = cond.getMach();
610 // At subsonic speeds CP at quarter chord
614 // At supersonic speeds use empirical formula
615 double beta = cond.getBeta();
616 return (ar * beta - 0.67) / (2 * ar * beta - 1);
619 // In between use interpolation polynomial
623 for (int i = 0; i < poly.length; i++) {
632 * Calculate CP position interpolation polynomial coefficients from the
633 * fin geometry. This is a fifth order polynomial that satisfies
642 * where f(M) = (ar*sqrt(M^2-1) - 0.67) / (2*ar*sqrt(M^2-1) - 1).
644 * The values were calculated analytically in Mathematica. The coefficients
645 * are used as poly[0] + poly[1]*x + poly[2]*x^2 + ...
647 private void calculatePoly() {
648 double denom = pow2(1 - 3.4641 * ar); // common denominator
650 poly[5] = (-1.58025 * (-0.728769 + ar) * (-0.192105 + ar)) / denom;
651 poly[4] = (12.8395 * (-0.725688 + ar) * (-0.19292 + ar)) / denom;
652 poly[3] = (-39.5062 * (-0.72074 + ar) * (-0.194245 + ar)) / denom;
653 poly[2] = (55.3086 * (-0.711482 + ar) * (-0.196772 + ar)) / denom;
654 poly[1] = (-31.6049 * (-0.705375 + ar) * (-0.198476 + ar)) / denom;
655 poly[0] = (9.16049 * (-0.588838 + ar) * (-0.20624 + ar)) / denom;
659 // @SuppressWarnings("null")
660 // public static void main(String arg[]) {
661 // Rocket rocket = TestRocket.makeRocket();
662 // FinSet finset = null;
664 // Iterator<RocketComponent> iter = rocket.deepIterator();
665 // while (iter.hasNext()) {
666 // RocketComponent c = iter.next();
667 // if (c instanceof FinSet) {
668 // finset = (FinSet)c;
673 // ((TrapezoidFinSet)finset).setHeight(0.10);
674 // ((TrapezoidFinSet)finset).setRootChord(0.10);
675 // ((TrapezoidFinSet)finset).setTipChord(0.10);
676 // ((TrapezoidFinSet)finset).setSweep(0.0);
679 // FinSetCalc calc = new FinSetCalc(finset);
681 // calc.calculateFinGeometry();
682 // FlightConditions cond = new FlightConditions(new Configuration(rocket));
683 // for (double m=0; m < 3; m+=0.05) {
685 // cond.setAOA(0.0*Math.PI/180);
686 // double cna = calc.calculateFinCNa1(cond);
687 // System.out.printf("%5.2f "+cna+"\n", m);
694 public double calculatePressureDragForce(FlightConditions conditions,
695 double stagnationCD, double baseCD, WarningSet warnings) {
697 double mach = conditions.getMach();
700 // Pressure fore-drag
701 if (crossSection == FinSet.CrossSection.AIRFOIL ||
702 crossSection == FinSet.CrossSection.ROUNDED) {
704 // Round leading edge
706 drag = Math.pow(1 - pow2(mach), -0.417) - 1;
707 } else if (mach < 1) {
708 drag = 1 - 1.785 * (mach - 0.9);
710 drag = 1.214 - 0.502 / pow2(mach) + 0.1095 / pow2(pow2(mach));
713 } else if (crossSection == FinSet.CrossSection.SQUARE) {
716 throw new UnsupportedOperationException("Unsupported fin profile: " + crossSection);
719 // Slanted leading edge
720 drag *= pow2(cosGammaLead);
722 // Trailing edge drag
723 if (crossSection == FinSet.CrossSection.SQUARE) {
725 } else if (crossSection == FinSet.CrossSection.ROUNDED) {
728 // Airfoil assumed to have zero base drag
731 // Scale to correct reference area
732 drag *= finCount * span * thickness / conditions.getRefArea();
738 private void calculateInterferenceFinCount(FinSet component) {
739 RocketComponent parent = component.getParent();
740 if (parent == null) {
741 throw new IllegalStateException("fin set without parent component");
744 double lead = component.toRelative(Coordinate.NUL, parent)[0].x;
745 double trail = component.toRelative(new Coordinate(component.getLength()),
749 * The counting fails if the fin root chord is very small, in that case assume
750 * no other fin interference than this fin set.
752 if (trail - lead < 0.007) {
753 interferenceFinCount = finCount;
755 interferenceFinCount = 0;
756 for (RocketComponent c : parent.getChildren()) {
757 if (c instanceof FinSet) {
758 double finLead = c.toRelative(Coordinate.NUL, parent)[0].x;
759 double finTrail = c.toRelative(new Coordinate(c.getLength()), parent)[0].x;
761 // Compute overlap of the fins
763 if ((finLead < trail - 0.005) && (finTrail > lead + 0.005)) {
764 interferenceFinCount += ((FinSet) c).getFinCount();
769 if (interferenceFinCount < component.getFinCount()) {
770 throw new BugException("Counted " + interferenceFinCount + " parallel fins, " +
771 "when component itself has " + component.getFinCount() +
772 ", fin points=" + Arrays.toString(component.getFinPoints()));