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 = calculateNonAxialForces(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 = calculateNonAxialForces(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 = calculateNonAxialForces(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());
151 * Perform the actual CP calculation.
153 private AerodynamicForces calculateNonAxialForces(Configuration configuration, FlightConditions conditions,
154 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
156 checkCache(configuration);
158 AerodynamicForces total = new AerodynamicForces();
161 double radius = 0; // aft radius of previous component
162 double componentX = 0; // aft coordinate of previous component
163 AerodynamicForces forces = new AerodynamicForces();
165 if (warnings == null)
166 warnings = ignoreWarningSet;
168 if (conditions.getAOA() > 17.5 * Math.PI / 180)
169 warnings.add(new Warning.LargeAOA(conditions.getAOA()));
173 buildCalcMap(configuration);
175 for (RocketComponent component : configuration) {
177 // Skip non-aerodynamic components
178 if (!component.isAerodynamic())
181 // Check for discontinuities
182 if (component instanceof SymmetricComponent) {
183 SymmetricComponent sym = (SymmetricComponent) component;
184 // TODO:LOW: Ignores other cluster components (not clusterable)
185 double x = component.toAbsolute(Coordinate.NUL)[0].x;
187 // Check for lengthwise discontinuity
188 if (x > componentX + 0.0001) {
189 if (!MathUtil.equals(radius, 0)) {
190 warnings.add(Warning.DISCONTINUITY);
194 componentX = component.toAbsolute(new Coordinate(component.getLength()))[0].x;
196 // Check for radius discontinuity
197 if (!MathUtil.equals(sym.getForeRadius(), radius)) {
198 warnings.add(Warning.DISCONTINUITY);
199 // TODO: MEDIUM: Apply correction to values to cp and to map
201 radius = sym.getAftRadius();
204 // Call calculation method
206 calcMap.get(component).calculateNonaxialForces(conditions, forces, warnings);
207 forces.setCP(component.toAbsolute(forces.getCP())[0]);
208 forces.setCm(forces.getCN() * forces.getCP().x / conditions.getRefLength());
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 // Calculate the skin friction coefficient (assume non-roughness limited)
259 if (configuration.getRocket().isPerfectFinish()) {
261 // Assume partial laminar layer. Roughness-limitation is checked later.
265 } else if (Re < 5.39e5) {
267 Cf = 1.328 / MathUtil.safeSqrt(Re);
270 Cf = 1.0 / pow2(1.50 * Math.log(Re) - 5.6) - 1700 / Re;
273 // Compressibility correction
276 // Below Re=1e6 no correction
279 c1 = 1 - 0.1 * pow2(mach) * (Re - 1e6) / 2e6; // transition to turbulent
281 c1 = 1 - 0.1 * pow2(mach);
288 c2 = 1 + (1.0 / Math.pow(1 + 0.045 * pow2(mach), 0.25) - 1) * (Re - 1e6) / 2e6;
290 c2 = 1.0 / Math.pow(1 + 0.045 * pow2(mach), 0.25);
295 // Applying continuously around Mach 1
298 } else if (mach < 1.1) {
299 Cf *= (c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2);
307 // Assume fully turbulent. Roughness-limitation is checked later.
313 Cf = 1.0 / pow2(1.50 * Math.log(Re) - 5.6);
316 // Compressibility correction
319 c1 = 1 - 0.1 * pow2(mach);
322 c2 = 1 / Math.pow(1 + 0.15 * pow2(mach), 0.58);
324 // Applying continuously around Mach 1
327 } else if (mach < 1.1) {
328 Cf *= c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2;
335 // Roughness-limited value correction term
336 double roughnessCorrection;
338 roughnessCorrection = 1 - 0.1 * pow2(mach);
339 } else if (mach > 1.1) {
340 roughnessCorrection = 1 / (1 + 0.18 * pow2(mach));
342 c1 = 1 - 0.1 * pow2(0.9);
343 c2 = 1.0 / (1 + 0.18 * pow2(1.1));
344 roughnessCorrection = c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2;
350 * Calculate the friction drag coefficient.
352 * The body wetted area is summed up and finally corrected with the rocket
353 * fineness ratio (calculated in the same iteration). The fins are corrected
354 * for thickness as we go on.
357 double finFriction = 0;
358 double bodyFriction = 0;
359 double maxR = 0, len = 0;
361 double[] roughnessLimited = new double[Finish.values().length];
362 Arrays.fill(roughnessLimited, Double.NaN);
364 for (RocketComponent c : configuration) {
366 // Consider only SymmetricComponents and FinSets:
367 if (!(c instanceof SymmetricComponent) &&
368 !(c instanceof FinSet))
371 // Calculate the roughness-limited friction coefficient
372 Finish finish = ((ExternalComponent) c).getFinish();
373 if (Double.isNaN(roughnessLimited[finish.ordinal()])) {
374 roughnessLimited[finish.ordinal()] =
375 0.032 * Math.pow(finish.getRoughnessSize() / configuration.getLength(), 0.2) *
380 * Actual Cf is maximum of Cf and the roughness-limited value.
381 * For perfect finish require additionally that Re > 1e6
384 if (configuration.getRocket().isPerfectFinish()) {
386 // For perfect finish require Re > 1e6
387 if ((Re > 1.0e6) && (roughnessLimited[finish.ordinal()] > Cf)) {
388 componentCf = roughnessLimited[finish.ordinal()];
395 // For fully turbulent use simple max
396 componentCf = Math.max(Cf, roughnessLimited[finish.ordinal()]);
402 // Calculate the friction drag:
403 if (c instanceof SymmetricComponent) {
405 SymmetricComponent s = (SymmetricComponent) c;
407 bodyFriction += componentCf * s.getComponentWetArea();
411 map.get(c).setFrictionCD(componentCf * s.getComponentWetArea()
412 / conditions.getRefArea());
415 double r = Math.max(s.getForeRadius(), s.getAftRadius());
418 len += c.getLength();
420 } else if (c instanceof FinSet) {
422 FinSet f = (FinSet) c;
423 double mac = ((FinSetCalc) calcMap.get(c)).getMACLength();
424 double cd = componentCf * (1 + 2 * f.getThickness() / mac) *
425 2 * f.getFinCount() * f.getFinArea();
429 map.get(c).setFrictionCD(cd / conditions.getRefArea());
435 // fB may be POSITIVE_INFINITY, but that's ok for us
436 double fB = (len + 0.0001) / maxR;
437 double correction = (1 + 1.0 / (2 * fB));
439 // Correct body data in map
441 for (RocketComponent c : map.keySet()) {
442 if (c instanceof SymmetricComponent) {
443 map.get(c).setFrictionCD(map.get(c).getFrictionCD() * correction);
448 return (finFriction + correction * bodyFriction) / conditions.getRefArea();
453 private double calculatePressureDrag(Configuration configuration, FlightConditions conditions,
454 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
456 double stagnation, base, total;
460 buildCalcMap(configuration);
462 stagnation = calculateStagnationCD(conditions.getMach());
463 base = calculateBaseCD(conditions.getMach());
466 for (RocketComponent c : configuration) {
467 if (!c.isAerodynamic())
470 // Pressure fore drag
471 double cd = calcMap.get(c).calculatePressureDragForce(conditions, stagnation, base,
476 map.get(c).setPressureCD(cd);
481 if (c instanceof SymmetricComponent) {
482 SymmetricComponent s = (SymmetricComponent) c;
484 if (radius < s.getForeRadius()) {
485 double area = Math.PI * (pow2(s.getForeRadius()) - pow2(radius));
486 cd = stagnation * area / conditions.getRefArea();
489 map.get(c).setPressureCD(map.get(c).getPressureCD() + cd);
493 radius = s.getAftRadius();
501 private double calculateBaseDrag(Configuration configuration, FlightConditions conditions,
502 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
506 RocketComponent prevComponent = null;
509 buildCalcMap(configuration);
511 base = calculateBaseCD(conditions.getMach());
514 for (RocketComponent c : configuration) {
515 if (!(c instanceof SymmetricComponent))
518 SymmetricComponent s = (SymmetricComponent) c;
520 if (radius > s.getForeRadius()) {
521 double area = Math.PI * (pow2(radius) - pow2(s.getForeRadius()));
522 double cd = base * area / conditions.getRefArea();
525 map.get(prevComponent).setBaseCD(cd);
529 radius = s.getAftRadius();
534 double area = Math.PI * pow2(radius);
535 double cd = base * area / conditions.getRefArea();
538 map.get(prevComponent).setBaseCD(cd);
547 public static double calculateStagnationCD(double m) {
550 pressure = 1 + pow2(m) / 4 + pow2(pow2(m)) / 40;
552 pressure = 1.84 - 0.76 / pow2(m) + 0.166 / pow2(pow2(m)) + 0.035 / pow2(m * m * m);
554 return 0.85 * pressure;
558 public static double calculateBaseCD(double m) {
560 return 0.12 + 0.13 * m * m;
568 private static final double[] axialDragPoly1, axialDragPoly2;
570 PolyInterpolator interpolator;
571 interpolator = new PolyInterpolator(
572 new double[] { 0, 17 * Math.PI / 180 },
573 new double[] { 0, 17 * Math.PI / 180 }
575 axialDragPoly1 = interpolator.interpolator(1, 1.3, 0, 0);
577 interpolator = new PolyInterpolator(
578 new double[] { 17 * Math.PI / 180, Math.PI / 2 },
579 new double[] { 17 * Math.PI / 180, Math.PI / 2 },
580 new double[] { Math.PI / 2 }
582 axialDragPoly2 = interpolator.interpolator(1.3, 0, 0, 0, 0);
587 * Calculate the axial drag from the total drag coefficient.
593 private double calculateAxialDrag(FlightConditions conditions, double cd) {
594 double aoa = MathUtil.clamp(conditions.getAOA(), 0, Math.PI);
597 // double sinaoa = conditions.getSinAOA();
598 // return cd * (1 + Math.min(sinaoa, 0.25));
601 if (aoa > Math.PI / 2)
603 if (aoa < 17 * Math.PI / 180)
604 mul = PolyInterpolator.eval(aoa, axialDragPoly1);
606 mul = PolyInterpolator.eval(aoa, axialDragPoly2);
608 if (conditions.getAOA() < Math.PI / 2)
615 private void calculateDampingMoments(Configuration configuration, FlightConditions conditions,
616 AerodynamicForces total) {
618 // Calculate pitch and yaw damping moments
619 double mul = getDampingMultiplier(configuration, conditions,
620 conditions.getPitchCenter().x);
621 double pitch = conditions.getPitchRate();
622 double yaw = conditions.getYawRate();
623 double vel = conditions.getVelocity();
625 vel = MathUtil.max(vel, 1);
627 mul *= 3; // TODO: Higher damping yields much more realistic apogee turn
629 total.setPitchDampingMoment(mul * MathUtil.sign(pitch) * pow2(pitch / vel));
630 total.setYawDampingMoment(mul * MathUtil.sign(yaw) * pow2(yaw / vel));
633 // TODO: MEDIUM: Are the rotation etc. being added correctly? sin/cos theta?
636 private double getDampingMultiplier(Configuration configuration, FlightConditions conditions,
638 if (cacheDiameter < 0) {
643 for (RocketComponent c : configuration) {
644 if (c instanceof SymmetricComponent) {
645 SymmetricComponent s = (SymmetricComponent) c;
646 area += s.getComponentPlanformArea();
647 cacheLength += s.getLength();
651 cacheDiameter = area / cacheLength;
657 mul = 0.275 * cacheDiameter / (conditions.getRefArea() * conditions.getRefLength());
658 mul *= (MathUtil.pow4(cgx) + MathUtil.pow4(cacheLength - cgx));
661 // TODO: LOW: This could be optimized a lot...
662 for (RocketComponent c : configuration) {
663 if (c instanceof FinSet) {
664 FinSet f = (FinSet) c;
665 mul += 0.6 * Math.min(f.getFinCount(), 4) * f.getFinArea() *
666 MathUtil.pow3(Math.abs(f.toAbsolute(new Coordinate(
667 ((FinSetCalc) calcMap.get(f)).getMidchordPos()))[0].x
669 (conditions.getRefArea() * conditions.getRefLength());
678 //////// The calculator map
681 protected void voidAerodynamicCache() {
682 super.voidAerodynamicCache();
690 private void buildCalcMap(Configuration configuration) {
691 Iterator<RocketComponent> iterator;
693 calcMap = new HashMap<RocketComponent, RocketComponentCalc>();
695 iterator = configuration.getRocket().iterator();
696 while (iterator.hasNext()) {
697 RocketComponent c = iterator.next();
699 if (!c.isAerodynamic())
702 calcMap.put(c, (RocketComponentCalc) Reflection.construct(BARROWMAN_PACKAGE,
703 c, BARROWMAN_SUFFIX, c));
709 public int getModID() {
710 // Only cached data is stored, return constant mod ID