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.ExternalComponent.Finish;
16 import net.sf.openrocket.rocketcomponent.FinSet;
17 import net.sf.openrocket.rocketcomponent.RocketComponent;
18 import net.sf.openrocket.rocketcomponent.SymmetricComponent;
19 import net.sf.openrocket.util.Coordinate;
20 import net.sf.openrocket.util.MathUtil;
21 import net.sf.openrocket.util.PolyInterpolator;
22 import net.sf.openrocket.util.Reflection;
25 * An aerodynamic calculator that uses the extended Barrowman method to
26 * calculate the CP of a rocket.
28 * @author Sampo Niskanen <sampo.niskanen@iki.fi>
30 public class BarrowmanCalculator extends AbstractAerodynamicCalculator {
32 private static final String BARROWMAN_PACKAGE = "net.sf.openrocket.aerodynamics.barrowman";
33 private static final String BARROWMAN_SUFFIX = "Calc";
36 private Map<RocketComponent, RocketComponentCalc> calcMap = null;
38 private double cacheDiameter = -1;
39 private double cacheLength = -1;
43 public BarrowmanCalculator() {
49 public BarrowmanCalculator newInstance() {
50 return new BarrowmanCalculator();
55 * Calculate the CP according to the extended Barrowman method.
58 public Coordinate getCP(Configuration configuration, FlightConditions conditions,
59 WarningSet warnings) {
60 checkCache(configuration);
61 AerodynamicForces forces = getNonAxial(configuration, conditions, null, warnings);
62 return forces.getCP();
68 public Map<RocketComponent, AerodynamicForces> getForceAnalysis(Configuration configuration,
69 FlightConditions conditions, WarningSet warnings) {
70 checkCache(configuration);
73 Map<RocketComponent, AerodynamicForces> map =
74 new LinkedHashMap<RocketComponent, AerodynamicForces>();
76 // Add all components to the map
77 for (RocketComponent c : configuration) {
78 f = new AerodynamicForces();
85 // Calculate non-axial force data
86 AerodynamicForces total = getNonAxial(configuration, conditions, map, warnings);
89 // Calculate friction data
90 total.setFrictionCD(calculateFrictionDrag(configuration, conditions, map, warnings));
91 total.setPressureCD(calculatePressureDrag(configuration, conditions, map, warnings));
92 total.setBaseCD(calculateBaseDrag(configuration, conditions, map, warnings));
94 total.setComponent(configuration.getRocket());
95 map.put(total.getComponent(), total);
98 for (RocketComponent c : map.keySet()) {
100 if (Double.isNaN(f.getBaseCD()) && Double.isNaN(f.getPressureCD()) &&
101 Double.isNaN(f.getFrictionCD()))
103 if (Double.isNaN(f.getBaseCD()))
105 if (Double.isNaN(f.getPressureCD()))
107 if (Double.isNaN(f.getFrictionCD()))
109 f.setCD(f.getBaseCD() + f.getPressureCD() + f.getFrictionCD());
110 f.setCaxial(calculateAxialDrag(conditions, f.getCD()));
119 public AerodynamicForces getAerodynamicForces(Configuration configuration,
120 FlightConditions conditions, WarningSet warnings) {
121 checkCache(configuration);
123 if (warnings == null)
124 warnings = ignoreWarningSet;
126 // Calculate non-axial force data
127 AerodynamicForces total = getNonAxial(configuration, conditions, null, warnings);
129 // Calculate friction data
130 total.setFrictionCD(calculateFrictionDrag(configuration, conditions, null, warnings));
131 total.setPressureCD(calculatePressureDrag(configuration, conditions, null, warnings));
132 total.setBaseCD(calculateBaseDrag(configuration, conditions, null, warnings));
134 total.setCD(total.getFrictionCD() + total.getPressureCD() + total.getBaseCD());
136 total.setCaxial(calculateAxialDrag(conditions, total.getCD()));
138 // Calculate pitch and yaw damping moments
139 calculateDampingMoments(configuration, conditions, total);
140 total.setCm(total.getCm() - total.getPitchDampingMoment());
141 total.setCyaw(total.getCyaw() - total.getYawDampingMoment());
150 * Perform the actual CP calculation.
152 private AerodynamicForces getNonAxial(Configuration configuration, FlightConditions conditions,
153 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
155 checkCache(configuration);
157 AerodynamicForces total = new AerodynamicForces();
160 double radius = 0; // aft radius of previous component
161 double componentX = 0; // aft coordinate of previous component
162 AerodynamicForces forces = new AerodynamicForces();
164 if (warnings == null)
165 warnings = ignoreWarningSet;
167 if (conditions.getAOA() > 17.5 * Math.PI / 180)
168 warnings.add(new Warning.LargeAOA(conditions.getAOA()));
172 buildCalcMap(configuration);
174 for (RocketComponent component : configuration) {
176 // Skip non-aerodynamic components
177 if (!component.isAerodynamic())
180 // Check for discontinuities
181 if (component instanceof SymmetricComponent) {
182 SymmetricComponent sym = (SymmetricComponent) component;
183 // TODO:LOW: Ignores other cluster components (not clusterable)
184 double x = component.toAbsolute(Coordinate.NUL)[0].x;
186 // Check for lengthwise discontinuity
187 if (x > componentX + 0.0001) {
188 if (!MathUtil.equals(radius, 0)) {
189 warnings.add(Warning.DISCONTINUITY);
193 componentX = component.toAbsolute(new Coordinate(component.getLength()))[0].x;
195 // Check for radius discontinuity
196 if (!MathUtil.equals(sym.getForeRadius(), radius)) {
197 warnings.add(Warning.DISCONTINUITY);
198 // TODO: MEDIUM: Apply correction to values to cp and to map
200 radius = sym.getAftRadius();
203 // Call calculation method
205 calcMap.get(component).calculateNonaxialForces(conditions, forces, warnings);
206 forces.setCP(component.toAbsolute(forces.getCP())[0]);
207 forces.setCm(forces.getCN() * forces.getCP().x / conditions.getRefLength());
208 // System.out.println(" CN="+forces.CN+" cp.x="+forces.cp.x+" Cm="+forces.Cm);
211 AerodynamicForces f = map.get(component);
213 f.setCP(forces.getCP());
214 f.setCNa(forces.getCNa());
215 f.setCN(forces.getCN());
216 f.setCm(forces.getCm());
217 f.setCside(forces.getCside());
218 f.setCyaw(forces.getCyaw());
219 f.setCroll(forces.getCroll());
220 f.setCrollDamp(forces.getCrollDamp());
221 f.setCrollForce(forces.getCrollForce());
224 total.setCP(total.getCP().average(forces.getCP()));
225 total.setCNa(total.getCNa() + forces.getCNa());
226 total.setCN(total.getCN() + forces.getCN());
227 total.setCm(total.getCm() + forces.getCm());
228 total.setCside(total.getCside() + forces.getCside());
229 total.setCyaw(total.getCyaw() + forces.getCyaw());
230 total.setCroll(total.getCroll() + forces.getCroll());
231 total.setCrollDamp(total.getCrollDamp() + forces.getCrollDamp());
232 total.setCrollForce(total.getCrollForce() + forces.getCrollForce());
241 //////////////// DRAG CALCULATIONS ////////////////
244 private double calculateFrictionDrag(Configuration configuration, FlightConditions conditions,
245 Map<RocketComponent, AerodynamicForces> map, WarningSet set) {
246 double c1 = 1.0, c2 = 1.0;
248 double mach = conditions.getMach();
253 buildCalcMap(configuration);
255 Re = conditions.getVelocity() * configuration.getLength() /
256 conditions.getAtmosphericConditions().getKinematicViscosity();
258 // System.out.printf("Re=%.3e ", Re);
260 // Calculate the skin friction coefficient (assume non-roughness limited)
261 if (configuration.getRocket().isPerfectFinish()) {
263 // System.out.printf("Perfect finish: Re=%f ",Re);
264 // Assume partial laminar layer. Roughness-limitation is checked later.
268 // System.out.printf("constant Cf=%f ",Cf);
269 } else if (Re < 5.39e5) {
271 Cf = 1.328 / MathUtil.safeSqrt(Re);
272 // System.out.printf("basic Cf=%f ",Cf);
275 Cf = 1.0 / pow2(1.50 * Math.log(Re) - 5.6) - 1700 / Re;
276 // System.out.printf("transitional Cf=%f ",Cf);
279 // Compressibility correction
282 // Below Re=1e6 no correction
285 c1 = 1 - 0.1 * pow2(mach) * (Re - 1e6) / 2e6; // transition to turbulent
287 c1 = 1 - 0.1 * pow2(mach);
294 c2 = 1 + (1.0 / Math.pow(1 + 0.045 * pow2(mach), 0.25) - 1) * (Re - 1e6) / 2e6;
296 c2 = 1.0 / Math.pow(1 + 0.045 * pow2(mach), 0.25);
301 // System.out.printf("c1=%f c2=%f\n", c1,c2);
302 // Applying continuously around Mach 1
305 } else if (mach < 1.1) {
306 Cf *= (c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2);
311 // System.out.printf("M=%f Cf=%f (smooth)\n",mach,Cf);
315 // Assume fully turbulent. Roughness-limitation is checked later.
319 // System.out.printf("LOW-TURB ");
322 Cf = 1.0 / pow2(1.50 * Math.log(Re) - 5.6);
323 // System.out.printf("NORMAL-TURB ");
326 // Compressibility correction
329 c1 = 1 - 0.1 * pow2(mach);
332 c2 = 1 / Math.pow(1 + 0.15 * pow2(mach), 0.58);
334 // Applying continuously around Mach 1
337 } else if (mach < 1.1) {
338 Cf *= c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2;
343 // System.out.printf("M=%f, Cd=%f (turbulent)\n", mach,Cf);
347 // Roughness-limited value correction term
348 double roughnessCorrection;
350 roughnessCorrection = 1 - 0.1 * pow2(mach);
351 } else if (mach > 1.1) {
352 roughnessCorrection = 1 / (1 + 0.18 * pow2(mach));
354 c1 = 1 - 0.1 * pow2(0.9);
355 c2 = 1.0 / (1 + 0.18 * pow2(1.1));
356 roughnessCorrection = c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2;
359 // System.out.printf("Cf=%.3f ", Cf);
363 * Calculate the friction drag coefficient.
365 * The body wetted area is summed up and finally corrected with the rocket
366 * fineness ratio (calculated in the same iteration). The fins are corrected
367 * for thickness as we go on.
370 double finFriction = 0;
371 double bodyFriction = 0;
372 double maxR = 0, len = 0;
374 double[] roughnessLimited = new double[Finish.values().length];
375 Arrays.fill(roughnessLimited, Double.NaN);
377 for (RocketComponent c : configuration) {
379 // Consider only SymmetricComponents and FinSets:
380 if (!(c instanceof SymmetricComponent) &&
381 !(c instanceof FinSet))
384 // Calculate the roughness-limited friction coefficient
385 Finish finish = ((ExternalComponent) c).getFinish();
386 if (Double.isNaN(roughnessLimited[finish.ordinal()])) {
387 roughnessLimited[finish.ordinal()] =
388 0.032 * Math.pow(finish.getRoughnessSize() / configuration.getLength(), 0.2) *
391 // System.out.printf("roughness["+finish+"]=%.3f ",
392 // roughnessLimited[finish.ordinal()]);
396 * Actual Cf is maximum of Cf and the roughness-limited value.
397 * For perfect finish require additionally that Re > 1e6
400 if (configuration.getRocket().isPerfectFinish()) {
402 // For perfect finish require Re > 1e6
403 if ((Re > 1.0e6) && (roughnessLimited[finish.ordinal()] > Cf)) {
404 componentCf = roughnessLimited[finish.ordinal()];
405 // System.out.printf(" rl=%f Cf=%f (perfect=%b)\n",
406 // roughnessLimited[finish.ordinal()],
407 // Cf,rocket.isPerfectFinish());
409 // System.out.printf("LIMITED ");
412 // System.out.printf("NORMAL ");
417 // For fully turbulent use simple max
418 componentCf = Math.max(Cf, roughnessLimited[finish.ordinal()]);
422 // System.out.printf("compCf=%.3f ", componentCf);
427 // Calculate the friction drag:
428 if (c instanceof SymmetricComponent) {
430 SymmetricComponent s = (SymmetricComponent) c;
432 bodyFriction += componentCf * s.getComponentWetArea();
436 map.get(c).setFrictionCD(componentCf * s.getComponentWetArea()
437 / conditions.getRefArea());
440 double r = Math.max(s.getForeRadius(), s.getAftRadius());
443 len += c.getLength();
445 } else if (c instanceof FinSet) {
447 FinSet f = (FinSet) c;
448 double mac = ((FinSetCalc) calcMap.get(c)).getMACLength();
449 double cd = componentCf * (1 + 2 * f.getThickness() / mac) *
450 2 * f.getFinCount() * f.getFinArea();
454 map.get(c).setFrictionCD(cd / conditions.getRefArea());
460 // fB may be POSITIVE_INFINITY, but that's ok for us
461 double fB = (len + 0.0001) / maxR;
462 double correction = (1 + 1.0 / (2 * fB));
464 // Correct body data in map
466 for (RocketComponent c : map.keySet()) {
467 if (c instanceof SymmetricComponent) {
468 map.get(c).setFrictionCD(map.get(c).getFrictionCD() * correction);
473 // System.out.printf("\n");
474 return (finFriction + correction * bodyFriction) / conditions.getRefArea();
479 private double calculatePressureDrag(Configuration configuration, FlightConditions conditions,
480 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
482 double stagnation, base, total;
486 buildCalcMap(configuration);
488 stagnation = calculateStagnationCD(conditions.getMach());
489 base = calculateBaseCD(conditions.getMach());
492 for (RocketComponent c : configuration) {
493 if (!c.isAerodynamic())
496 // Pressure fore drag
497 double cd = calcMap.get(c).calculatePressureDragForce(conditions, stagnation, base,
502 map.get(c).setPressureCD(cd);
507 if (c instanceof SymmetricComponent) {
508 SymmetricComponent s = (SymmetricComponent) c;
510 if (radius < s.getForeRadius()) {
511 double area = Math.PI * (pow2(s.getForeRadius()) - pow2(radius));
512 cd = stagnation * area / conditions.getRefArea();
515 map.get(c).setPressureCD(map.get(c).getPressureCD() + cd);
519 radius = s.getAftRadius();
527 private double calculateBaseDrag(Configuration configuration, FlightConditions conditions,
528 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
532 RocketComponent prevComponent = null;
535 buildCalcMap(configuration);
537 base = calculateBaseCD(conditions.getMach());
540 for (RocketComponent c : configuration) {
541 if (!(c instanceof SymmetricComponent))
544 SymmetricComponent s = (SymmetricComponent) c;
546 if (radius > s.getForeRadius()) {
547 double area = Math.PI * (pow2(radius) - pow2(s.getForeRadius()));
548 double cd = base * area / conditions.getRefArea();
551 map.get(prevComponent).setBaseCD(cd);
555 radius = s.getAftRadius();
560 double area = Math.PI * pow2(radius);
561 double cd = base * area / conditions.getRefArea();
564 map.get(prevComponent).setBaseCD(cd);
573 public static double calculateStagnationCD(double m) {
576 pressure = 1 + pow2(m) / 4 + pow2(pow2(m)) / 40;
578 pressure = 1.84 - 0.76 / pow2(m) + 0.166 / pow2(pow2(m)) + 0.035 / pow2(m * m * m);
580 return 0.85 * pressure;
584 public static double calculateBaseCD(double m) {
586 return 0.12 + 0.13 * m * m;
594 private static final double[] axialDragPoly1, axialDragPoly2;
596 PolyInterpolator interpolator;
597 interpolator = new PolyInterpolator(
598 new double[] { 0, 17 * Math.PI / 180 },
599 new double[] { 0, 17 * Math.PI / 180 }
601 axialDragPoly1 = interpolator.interpolator(1, 1.3, 0, 0);
603 interpolator = new PolyInterpolator(
604 new double[] { 17 * Math.PI / 180, Math.PI / 2 },
605 new double[] { 17 * Math.PI / 180, Math.PI / 2 },
606 new double[] { Math.PI / 2 }
608 axialDragPoly2 = interpolator.interpolator(1.3, 0, 0, 0, 0);
613 * Calculate the axial drag from the total drag coefficient.
619 private double calculateAxialDrag(FlightConditions conditions, double cd) {
620 double aoa = MathUtil.clamp(conditions.getAOA(), 0, Math.PI);
623 // double sinaoa = conditions.getSinAOA();
624 // return cd * (1 + Math.min(sinaoa, 0.25));
627 if (aoa > Math.PI / 2)
629 if (aoa < 17 * Math.PI / 180)
630 mul = PolyInterpolator.eval(aoa, axialDragPoly1);
632 mul = PolyInterpolator.eval(aoa, axialDragPoly2);
634 if (conditions.getAOA() < Math.PI / 2)
641 private void calculateDampingMoments(Configuration configuration, FlightConditions conditions,
642 AerodynamicForces total) {
644 // Calculate pitch and yaw damping moments
645 double mul = getDampingMultiplier(configuration, conditions,
646 conditions.getPitchCenter().x);
647 double pitch = conditions.getPitchRate();
648 double yaw = conditions.getYawRate();
649 double vel = conditions.getVelocity();
651 // double Cm = total.Cm - total.CN * total.cg.x / conditions.getRefLength();
652 // System.out.printf("Damping pitch/yaw, mul=%.4f pitch rate=%.4f "+
653 // "Cm=%.4f / %.4f effect=%.4f aoa=%.4f\n", mul, pitch, total.Cm, Cm,
654 // -(mul * MathUtil.sign(pitch) * pow2(pitch/vel)),
655 // conditions.getAOA()*180/Math.PI);
657 mul *= 3; // TODO: Higher damping yields much more realistic apogee turn
659 // total.Cm -= mul * pitch / pow2(vel);
660 // total.Cyaw -= mul * yaw / pow2(vel);
661 total.setPitchDampingMoment(mul * MathUtil.sign(pitch) * pow2(pitch / vel));
662 total.setYawDampingMoment(mul * MathUtil.sign(yaw) * pow2(yaw / vel));
666 // TODO: MEDIUM: Are the rotation etc. being added correctly? sin/cos theta?
669 private double getDampingMultiplier(Configuration configuration, FlightConditions conditions,
671 if (cacheDiameter < 0) {
676 for (RocketComponent c : configuration) {
677 if (c instanceof SymmetricComponent) {
678 SymmetricComponent s = (SymmetricComponent) c;
679 area += s.getComponentPlanformArea();
680 cacheLength += s.getLength();
684 cacheDiameter = area / cacheLength;
690 mul = 0.275 * cacheDiameter / (conditions.getRefArea() * conditions.getRefLength());
691 mul *= (MathUtil.pow4(cgx) + MathUtil.pow4(cacheLength - cgx));
694 // TODO: LOW: This could be optimized a lot...
695 for (RocketComponent c : configuration) {
696 if (c instanceof FinSet) {
697 FinSet f = (FinSet) c;
698 mul += 0.6 * Math.min(f.getFinCount(), 4) * f.getFinArea() *
699 MathUtil.pow3(Math.abs(f.toAbsolute(new Coordinate(
700 ((FinSetCalc) calcMap.get(f)).getMidchordPos()))[0].x
702 (conditions.getRefArea() * conditions.getRefLength());
711 //////// The calculator map
714 protected void voidAerodynamicCache() {
715 super.voidAerodynamicCache();
723 private void buildCalcMap(Configuration configuration) {
724 Iterator<RocketComponent> iterator;
726 calcMap = new HashMap<RocketComponent, RocketComponentCalc>();
728 iterator = configuration.getRocket().iterator();
729 while (iterator.hasNext()) {
730 RocketComponent c = iterator.next();
732 if (!c.isAerodynamic())
735 calcMap.put(c, (RocketComponentCalc) Reflection.construct(BARROWMAN_PACKAGE,
736 c, BARROWMAN_SUFFIX, c));
742 public int getModID() {
743 // Only cached data is stored, return constant mod ID