1 package net.sf.openrocket.aerodynamics;
3 import static net.sf.openrocket.util.MathUtil.pow2;
5 import java.util.Arrays;
6 import java.util.HashMap;
7 import java.util.Iterator;
8 import java.util.LinkedHashMap;
11 import net.sf.openrocket.aerodynamics.barrowman.FinSetCalc;
12 import net.sf.openrocket.aerodynamics.barrowman.RocketComponentCalc;
13 import net.sf.openrocket.rocketcomponent.Configuration;
14 import net.sf.openrocket.rocketcomponent.ExternalComponent;
15 import net.sf.openrocket.rocketcomponent.FinSet;
16 import net.sf.openrocket.rocketcomponent.Rocket;
17 import net.sf.openrocket.rocketcomponent.RocketComponent;
18 import net.sf.openrocket.rocketcomponent.SymmetricComponent;
19 import net.sf.openrocket.rocketcomponent.ExternalComponent.Finish;
20 import net.sf.openrocket.util.Coordinate;
21 import net.sf.openrocket.util.MathUtil;
22 import net.sf.openrocket.util.PolyInterpolator;
23 import net.sf.openrocket.util.Reflection;
24 import net.sf.openrocket.util.Test;
27 * An aerodynamic calculator that uses the extended Barrowman method to
28 * calculate the CP of a rocket.
30 * @author Sampo Niskanen <sampo.niskanen@iki.fi>
32 public class BarrowmanCalculator extends AerodynamicCalculator {
34 private static final String BARROWMAN_PACKAGE = "net.sf.openrocket.aerodynamics.barrowman";
35 private static final String BARROWMAN_SUFFIX = "Calc";
38 private Map<RocketComponent, RocketComponentCalc> calcMap = null;
40 private double cacheDiameter = -1;
41 private double cacheLength = -1;
45 public BarrowmanCalculator() {
49 public BarrowmanCalculator(Configuration config) {
55 public BarrowmanCalculator newInstance() {
56 return new BarrowmanCalculator();
61 * Calculate the CP according to the extended Barrowman method.
64 public Coordinate getCP(FlightConditions conditions, WarningSet warnings) {
65 AerodynamicForces forces = getNonAxial(conditions, null, warnings);
72 public Map<RocketComponent, AerodynamicForces> getForceAnalysis(FlightConditions conditions,
73 WarningSet warnings) {
76 Map<RocketComponent, AerodynamicForces> map =
77 new LinkedHashMap<RocketComponent, AerodynamicForces>();
79 // Add all components to the map
80 for (RocketComponent c: configuration) {
81 f = new AerodynamicForces();
85 f.cg = Coordinate.NUL;
86 for (Coordinate coord: c.toAbsolute(c.getCG())) {
87 f.cg = f.cg.average(coord);
94 // Calculate non-axial force data
95 AerodynamicForces total = getNonAxial(conditions, map, warnings);
98 // Calculate friction data
99 total.frictionCD = calculateFrictionDrag(conditions, map, warnings);
100 total.pressureCD = calculatePressureDrag(conditions, map, warnings);
101 total.baseCD = calculateBaseDrag(conditions, map, warnings);
104 total.component = rocket;
105 map.put(rocket, total);
108 for (RocketComponent c: map.keySet()) {
110 if (Double.isNaN(f.baseCD) && Double.isNaN(f.pressureCD) &&
111 Double.isNaN(f.frictionCD))
113 if (Double.isNaN(f.baseCD))
115 if (Double.isNaN(f.pressureCD))
117 if (Double.isNaN(f.frictionCD))
119 f.CD = f.baseCD + f.pressureCD + f.frictionCD;
120 f.Caxial = calculateAxialDrag(conditions, f.CD);
129 public AerodynamicForces getAerodynamicForces(double time, FlightConditions conditions,
130 WarningSet warnings) {
132 if (warnings == null)
133 warnings = ignoreWarningSet;
135 // Calculate non-axial force data
136 AerodynamicForces total = getNonAxial(conditions, null, warnings);
138 // Calculate friction data
139 total.frictionCD = calculateFrictionDrag(conditions, null, warnings);
140 total.pressureCD = calculatePressureDrag(conditions, null, warnings);
141 total.baseCD = calculateBaseDrag(conditions, null, warnings);
143 total.CD = total.frictionCD + total.pressureCD + total.baseCD;
145 total.Caxial = calculateAxialDrag(conditions, total.CD);
148 // Calculate CG and moments of inertia
149 total.cg = this.getCG(time);
150 total.longitudalInertia = this.getLongitudalInertia(time);
151 total.rotationalInertia = this.getRotationalInertia(time);
154 // Calculate pitch and yaw damping moments
155 calculateDampingMoments(conditions, total);
156 total.Cm -= total.pitchDampingMoment;
157 total.Cyaw -= total.yawDampingMoment;
160 // System.out.println("Conditions are "+conditions + "
161 // pitch rate="+conditions.getPitchRate());
162 // System.out.println("Total Cm="+total.Cm+" damping effect="+
163 // (12 * Math.signum(conditions.getPitchRate()) *
164 // MathUtil.pow2(conditions.getPitchRate()) /
165 // MathUtil.pow2(conditions.getVelocity())));
167 // double ef = Math.abs(12 *
168 // MathUtil.pow2(conditions.getPitchRate()) /
169 // MathUtil.pow2(conditions.getVelocity()));
171 //// System.out.println("maxEffect="+maxEffect);
172 // total.Cm -= 12 * Math.signum(conditions.getPitchRate()) *
173 // MathUtil.pow2(conditions.getPitchRate()) /
174 // MathUtil.pow2(conditions.getVelocity());
176 // total.Cyaw -= 0.06 * Math.signum(conditions.getYawRate()) *
177 // MathUtil.pow2(conditions.getYawRate()) /
178 // MathUtil.pow2(conditions.getVelocity());
186 public AerodynamicForces getAxialForces(double time,
187 FlightConditions conditions, WarningSet warnings) {
189 if (warnings == null)
190 warnings = ignoreWarningSet;
192 AerodynamicForces total = new AerodynamicForces();
195 // Calculate friction data
196 total.frictionCD = calculateFrictionDrag(conditions, null, warnings);
197 total.pressureCD = calculatePressureDrag(conditions, null, warnings);
198 total.baseCD = calculateBaseDrag(conditions, null, warnings);
200 total.CD = total.frictionCD + total.pressureCD + total.baseCD;
202 total.Caxial = calculateAxialDrag(conditions, total.CD);
204 // Calculate CG and moments of inertia
205 total.cg = this.getCG(time);
206 total.longitudalInertia = this.getLongitudalInertia(time);
207 total.rotationalInertia = this.getRotationalInertia(time);
217 * Perform the actual CP calculation.
219 private AerodynamicForces getNonAxial(FlightConditions conditions,
220 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
222 AerodynamicForces total = new AerodynamicForces();
225 double radius = 0; // aft radius of previous component
226 double componentX = 0; // aft coordinate of previous component
227 AerodynamicForces forces = new AerodynamicForces();
229 if (warnings == null)
230 warnings = ignoreWarningSet;
232 if (conditions.getAOA() > 17.5*Math.PI/180)
233 warnings.add(new Warning.LargeAOA(conditions.getAOA()));
240 for (RocketComponent component: configuration) {
242 // Skip non-aerodynamic components
243 if (!component.isAerodynamic())
246 // Check for discontinuities
247 if (component instanceof SymmetricComponent) {
248 SymmetricComponent sym = (SymmetricComponent) component;
249 // TODO:LOW: Ignores other cluster components (not clusterable)
250 double x = component.toAbsolute(Coordinate.NUL)[0].x;
252 // Check for lengthwise discontinuity
253 if (x > componentX + 0.0001){
254 if (!MathUtil.equals(radius, 0)) {
255 warnings.add(Warning.DISCONTINUITY);
259 componentX = component.toAbsolute(new Coordinate(component.getLength()))[0].x;
261 // Check for radius discontinuity
262 if (!MathUtil.equals(sym.getForeRadius(), radius)) {
263 warnings.add(Warning.DISCONTINUITY);
264 // TODO: MEDIUM: Apply correction to values to cp and to map
266 radius = sym.getAftRadius();
269 // Call calculation method
271 calcMap.get(component).calculateNonaxialForces(conditions, forces, warnings);
272 forces.cp = component.toAbsolute(forces.cp)[0];
273 forces.Cm = forces.CN * forces.cp.x / conditions.getRefLength();
274 // System.out.println(" CN="+forces.CN+" cp.x="+forces.cp.x+" Cm="+forces.Cm);
277 AerodynamicForces f = map.get(component);
283 f.Cside = forces.Cside;
284 f.Cyaw = forces.Cyaw;
285 f.Croll = forces.Croll;
286 f.CrollDamp = forces.CrollDamp;
287 f.CrollForce = forces.CrollForce;
290 total.cp = total.cp.average(forces.cp);
291 total.CNa += forces.CNa;
292 total.CN += forces.CN;
293 total.Cm += forces.Cm;
294 total.Cside += forces.Cside;
295 total.Cyaw += forces.Cyaw;
296 total.Croll += forces.Croll;
297 total.CrollDamp += forces.CrollDamp;
298 total.CrollForce += forces.CrollForce;
307 //////////////// DRAG CALCULATIONS ////////////////
310 private double calculateFrictionDrag(FlightConditions conditions,
311 Map<RocketComponent, AerodynamicForces> map, WarningSet set) {
312 double c1=1.0, c2=1.0;
314 double mach = conditions.getMach();
321 Re = conditions.getVelocity() * configuration.getLength() /
322 conditions.getAtmosphericConditions().getKinematicViscosity();
324 // System.out.printf("Re=%.3e ", Re);
326 // Calculate the skin friction coefficient (assume non-roughness limited)
327 if (configuration.getRocket().isPerfectFinish()) {
329 // System.out.printf("Perfect finish: Re=%f ",Re);
330 // Assume partial laminar layer. Roughness-limitation is checked later.
334 // System.out.printf("constant Cf=%f ",Cf);
335 } else if (Re < 5.39e5) {
337 Cf = 1.328 / Math.sqrt(Re);
338 // System.out.printf("basic Cf=%f ",Cf);
341 Cf = 1.0/pow2(1.50 * Math.log(Re) - 5.6) - 1700/Re;
342 // System.out.printf("transitional Cf=%f ",Cf);
345 // Compressibility correction
348 // Below Re=1e6 no correction
351 c1 = 1 - 0.1*pow2(mach)*(Re-1e6)/2e6; // transition to turbulent
353 c1 = 1 - 0.1*pow2(mach);
360 c2 = 1 + (1.0 / Math.pow(1+0.045*pow2(mach), 0.25) -1) * (Re-1e6)/2e6;
362 c2 = 1.0 / Math.pow(1+0.045*pow2(mach), 0.25);
367 // System.out.printf("c1=%f c2=%f\n", c1,c2);
368 // Applying continuously around Mach 1
371 } else if (mach < 1.1) {
372 Cf *= (c2 * (mach-0.9)/0.2 + c1 * (1.1-mach)/0.2);
377 // System.out.printf("M=%f Cf=%f (smooth)\n",mach,Cf);
381 // Assume fully turbulent. Roughness-limitation is checked later.
385 // System.out.printf("LOW-TURB ");
388 Cf = 1.0/pow2(1.50 * Math.log(Re) - 5.6);
389 // System.out.printf("NORMAL-TURB ");
392 // Compressibility correction
395 c1 = 1 - 0.1*pow2(mach);
398 c2 = 1/Math.pow(1 + 0.15*pow2(mach), 0.58);
400 // Applying continuously around Mach 1
403 } else if (mach < 1.1) {
404 Cf *= c2 * (mach-0.9)/0.2 + c1 * (1.1-mach)/0.2;
409 // System.out.printf("M=%f, Cd=%f (turbulent)\n", mach,Cf);
413 // Roughness-limited value correction term
414 double roughnessCorrection;
416 roughnessCorrection = 1 - 0.1*pow2(mach);
417 } else if (mach > 1.1) {
418 roughnessCorrection = 1/(1 + 0.18*pow2(mach));
420 c1 = 1 - 0.1*pow2(0.9);
421 c2 = 1.0/(1+0.18 * pow2(1.1));
422 roughnessCorrection = c2 * (mach-0.9)/0.2 + c1 * (1.1-mach)/0.2;
425 // System.out.printf("Cf=%.3f ", Cf);
429 * Calculate the friction drag coefficient.
431 * The body wetted area is summed up and finally corrected with the rocket
432 * fineness ratio (calculated in the same iteration). The fins are corrected
433 * for thickness as we go on.
436 double finFriction = 0;
437 double bodyFriction = 0;
438 double maxR=0, len=0;
440 double[] roughnessLimited = new double[Finish.values().length];
441 Arrays.fill(roughnessLimited, Double.NaN);
443 for (RocketComponent c: configuration) {
445 // Consider only SymmetricComponents and FinSets:
446 if (!(c instanceof SymmetricComponent) &&
447 !(c instanceof FinSet))
450 // Calculate the roughness-limited friction coefficient
451 Finish finish = ((ExternalComponent)c).getFinish();
452 if (Double.isNaN(roughnessLimited[finish.ordinal()])) {
453 roughnessLimited[finish.ordinal()] =
454 0.032 * Math.pow(finish.getRoughnessSize()/configuration.getLength(), 0.2) *
457 // System.out.printf("roughness["+finish+"]=%.3f ",
458 // roughnessLimited[finish.ordinal()]);
462 * Actual Cf is maximum of Cf and the roughness-limited value.
463 * For perfect finish require additionally that Re > 1e6
466 if (configuration.getRocket().isPerfectFinish()) {
468 // For perfect finish require Re > 1e6
469 if ((Re > 1.0e6) && (roughnessLimited[finish.ordinal()] > Cf)) {
470 componentCf = roughnessLimited[finish.ordinal()];
471 // System.out.printf(" rl=%f Cf=%f (perfect=%b)\n",
472 // roughnessLimited[finish.ordinal()],
473 // Cf,rocket.isPerfectFinish());
475 // System.out.printf("LIMITED ");
478 // System.out.printf("NORMAL ");
483 // For fully turbulent use simple max
484 componentCf = Math.max(Cf, roughnessLimited[finish.ordinal()]);
488 // System.out.printf("compCf=%.3f ", componentCf);
493 // Calculate the friction drag:
494 if (c instanceof SymmetricComponent) {
496 SymmetricComponent s = (SymmetricComponent)c;
498 bodyFriction += componentCf * s.getComponentWetArea();
502 map.get(c).frictionCD = componentCf * s.getComponentWetArea()
503 / conditions.getRefArea();
506 double r = Math.max(s.getForeRadius(), s.getAftRadius());
509 len += c.getLength();
511 } else if (c instanceof FinSet) {
513 FinSet f = (FinSet)c;
514 double mac = ((FinSetCalc)calcMap.get(c)).getMACLength();
515 double cd = componentCf * (1 + 2*f.getThickness()/mac) *
516 2*f.getFinCount() * f.getFinArea();
520 map.get(c).frictionCD = cd / conditions.getRefArea();
526 // fB may be POSITIVE_INFINITY, but that's ok for us
527 double fB = (len+0.0001) / maxR;
528 double correction = (1 + 1.0/(2*fB));
530 // Correct body data in map
532 for (RocketComponent c: map.keySet()) {
533 if (c instanceof SymmetricComponent) {
534 map.get(c).frictionCD *= correction;
539 // System.out.printf("\n");
540 return (finFriction + correction*bodyFriction) / conditions.getRefArea();
545 private double calculatePressureDrag(FlightConditions conditions,
546 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
548 double stagnation, base, total;
554 stagnation = calculateStagnationCD(conditions.getMach());
555 base = calculateBaseCD(conditions.getMach());
558 for (RocketComponent c: configuration) {
559 if (!c.isAerodynamic())
562 // Pressure fore drag
563 double cd = calcMap.get(c).calculatePressureDragForce(conditions, stagnation, base,
568 map.get(c).pressureCD = cd;
573 if (c instanceof SymmetricComponent) {
574 SymmetricComponent s = (SymmetricComponent)c;
576 if (radius < s.getForeRadius()) {
577 double area = Math.PI*(pow2(s.getForeRadius()) - pow2(radius));
578 cd = stagnation * area / conditions.getRefArea();
581 map.get(c).pressureCD += cd;
585 radius = s.getAftRadius();
593 private double calculateBaseDrag(FlightConditions conditions,
594 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
598 RocketComponent prevComponent = null;
603 base = calculateBaseCD(conditions.getMach());
606 for (RocketComponent c: configuration) {
607 if (!(c instanceof SymmetricComponent))
610 SymmetricComponent s = (SymmetricComponent)c;
612 if (radius > s.getForeRadius()) {
613 double area = Math.PI*(pow2(radius) - pow2(s.getForeRadius()));
614 double cd = base * area / conditions.getRefArea();
617 map.get(prevComponent).baseCD = cd;
621 radius = s.getAftRadius();
626 double area = Math.PI*pow2(radius);
627 double cd = base * area / conditions.getRefArea();
630 map.get(prevComponent).baseCD = cd;
639 public static double calculateStagnationCD(double m) {
642 pressure = 1 + pow2(m)/4 + pow2(pow2(m))/40;
644 pressure = 1.84 - 0.76/pow2(m) + 0.166/pow2(pow2(m)) + 0.035/pow2(m*m*m);
646 return 0.85 * pressure;
650 public static double calculateBaseCD(double m) {
652 return 0.12 + 0.13 * m*m;
660 private static final double[] axialDragPoly1, axialDragPoly2;
662 PolyInterpolator interpolator;
663 interpolator = new PolyInterpolator(
664 new double[] { 0, 17*Math.PI/180 },
665 new double[] { 0, 17*Math.PI/180 }
667 axialDragPoly1 = interpolator.interpolator(1, 1.3, 0, 0);
669 interpolator = new PolyInterpolator(
670 new double[] { 17*Math.PI/180, Math.PI/2 },
671 new double[] { 17*Math.PI/180, Math.PI/2 },
672 new double[] { Math.PI/2 }
674 axialDragPoly2 = interpolator.interpolator(1.3, 0, 0, 0, 0);
679 * Calculate the axial drag from the total drag coefficient.
685 private double calculateAxialDrag(FlightConditions conditions, double cd) {
686 double aoa = MathUtil.clamp(conditions.getAOA(), 0, Math.PI);
689 // double sinaoa = conditions.getSinAOA();
690 // return cd * (1 + Math.min(sinaoa, 0.25));
695 if (aoa < 17*Math.PI/180)
696 mul = PolyInterpolator.eval(aoa, axialDragPoly1);
698 mul = PolyInterpolator.eval(aoa, axialDragPoly2);
700 if (conditions.getAOA() < Math.PI/2)
707 private void calculateDampingMoments(FlightConditions conditions,
708 AerodynamicForces total) {
710 // Calculate pitch and yaw damping moments
711 if (conditions.getPitchRate() > 0.1 || conditions.getYawRate() > 0.1 || true) {
712 double mul = getDampingMultiplier(conditions, total.cg.x);
713 double pitch = conditions.getPitchRate();
714 double yaw = conditions.getYawRate();
715 double vel = conditions.getVelocity();
717 // double Cm = total.Cm - total.CN * total.cg.x / conditions.getRefLength();
718 // System.out.printf("Damping pitch/yaw, mul=%.4f pitch rate=%.4f "+
719 // "Cm=%.4f / %.4f effect=%.4f aoa=%.4f\n", mul, pitch, total.Cm, Cm,
720 // -(mul * MathUtil.sign(pitch) * pow2(pitch/vel)),
721 // conditions.getAOA()*180/Math.PI);
723 mul *= 3; // TODO: Higher damping yields much more realistic apogee turn
725 // total.Cm -= mul * pitch / pow2(vel);
726 // total.Cyaw -= mul * yaw / pow2(vel);
727 total.pitchDampingMoment = mul * MathUtil.sign(pitch) * pow2(pitch/vel);
728 total.yawDampingMoment = mul * MathUtil.sign(yaw) * pow2(yaw/vel);
730 total.pitchDampingMoment = 0;
731 total.yawDampingMoment = 0;
736 // TODO: MEDIUM: Are the rotation etc. being added correctly? sin/cos theta?
739 private double getDampingMultiplier(FlightConditions conditions, double cgx) {
740 if (cacheDiameter < 0) {
745 for (RocketComponent c: configuration) {
746 if (c instanceof SymmetricComponent) {
747 SymmetricComponent s = (SymmetricComponent)c;
748 area += s.getComponentPlanformArea();
749 cacheLength += s.getLength();
753 cacheDiameter = area / cacheLength;
759 mul = 0.275 * cacheDiameter / (conditions.getRefArea() * conditions.getRefLength());
760 mul *= (MathUtil.pow4(cgx) + MathUtil.pow4(cacheLength - cgx));
763 // TODO: LOW: This could be optimized a lot...
764 for (RocketComponent c: configuration) {
765 if (c instanceof FinSet) {
766 FinSet f = (FinSet)c;
767 mul += 0.6 * Math.min(f.getFinCount(), 4) * f.getFinArea() *
768 MathUtil.pow3(Math.abs(f.toAbsolute(new Coordinate(
769 ((FinSetCalc)calcMap.get(f)).getMidchordPos()))[0].x
771 (conditions.getRefArea() * conditions.getRefLength());
780 //////// The calculator map
783 protected void voidAerodynamicCache() {
784 super.voidAerodynamicCache();
792 private void buildCalcMap() {
793 Iterator<RocketComponent> iterator;
795 calcMap = new HashMap<RocketComponent, RocketComponentCalc>();
797 iterator = rocket.deepIterator();
798 while (iterator.hasNext()) {
799 RocketComponent c = iterator.next();
801 if (!c.isAerodynamic())
804 calcMap.put(c, (RocketComponentCalc) Reflection.construct(BARROWMAN_PACKAGE,
805 c, BARROWMAN_SUFFIX, c));
812 public static void main(String[] arg) {
814 PolyInterpolator interpolator;
816 interpolator = new PolyInterpolator(
817 new double[] { 0, 17*Math.PI/180 },
818 new double[] { 0, 17*Math.PI/180 }
820 double[] poly1 = interpolator.interpolator(1, 1.3, 0, 0);
822 interpolator = new PolyInterpolator(
823 new double[] { 17*Math.PI/180, Math.PI/2 },
824 new double[] { 17*Math.PI/180, Math.PI/2 },
825 new double[] { Math.PI/2 }
827 double[] poly2 = interpolator.interpolator(1.3, 0, 0, 0, 0);
830 for (double a=0; a<=180.1; a++) {
831 double r = a*Math.PI/180;
836 if (r < 18*Math.PI/180)
837 value = PolyInterpolator.eval(r, poly1);
839 value = PolyInterpolator.eval(r, poly2);
841 System.out.println(""+a+" "+value);
847 Rocket normal = Test.makeRocket();
848 Rocket perfect = Test.makeRocket();
849 normal.setPerfectFinish(false);
850 perfect.setPerfectFinish(true);
852 Configuration confNormal = new Configuration(normal);
853 Configuration confPerfect = new Configuration(perfect);
855 for (RocketComponent c: confNormal) {
856 if (c instanceof ExternalComponent) {
857 ((ExternalComponent)c).setFinish(Finish.NORMAL);
860 for (RocketComponent c: confPerfect) {
861 if (c instanceof ExternalComponent) {
862 ((ExternalComponent)c).setFinish(Finish.NORMAL);
867 confNormal.setToStage(0);
868 confPerfect.setToStage(0);
872 BarrowmanCalculator calcNormal = new BarrowmanCalculator(confNormal);
873 BarrowmanCalculator calcPerfect = new BarrowmanCalculator(confPerfect);
875 FlightConditions conditions = new FlightConditions(confNormal);
877 for (double mach=0; mach < 3; mach += 0.1) {
878 conditions.setMach(mach);
880 Map<RocketComponent, AerodynamicForces> data =
881 calcNormal.getForceAnalysis(conditions, null);
882 AerodynamicForces forcesNormal = data.get(normal);
884 data = calcPerfect.getForceAnalysis(conditions, null);
885 AerodynamicForces forcesPerfect = data.get(perfect);
887 System.out.printf("%f %f %f %f %f %f %f\n",mach,
888 forcesNormal.pressureCD, forcesPerfect.pressureCD,
889 forcesNormal.frictionCD, forcesPerfect.frictionCD,
890 forcesNormal.CD, forcesPerfect.CD);