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());
151 * Perform the actual CP calculation.
153 private AerodynamicForces getNonAxial(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());
209 // System.out.println(" CN="+forces.CN+" cp.x="+forces.cp.x+" Cm="+forces.Cm);
212 AerodynamicForces f = map.get(component);
214 f.setCP(forces.getCP());
215 f.setCNa(forces.getCNa());
216 f.setCN(forces.getCN());
217 f.setCm(forces.getCm());
218 f.setCside(forces.getCside());
219 f.setCyaw(forces.getCyaw());
220 f.setCroll(forces.getCroll());
221 f.setCrollDamp(forces.getCrollDamp());
222 f.setCrollForce(forces.getCrollForce());
225 total.setCP(total.getCP().average(forces.getCP()));
226 total.setCNa(total.getCNa() + forces.getCNa());
227 total.setCN(total.getCN() + forces.getCN());
228 total.setCm(total.getCm() + forces.getCm());
229 total.setCside(total.getCside() + forces.getCside());
230 total.setCyaw(total.getCyaw() + forces.getCyaw());
231 total.setCroll(total.getCroll() + forces.getCroll());
232 total.setCrollDamp(total.getCrollDamp() + forces.getCrollDamp());
233 total.setCrollForce(total.getCrollForce() + forces.getCrollForce());
242 //////////////// DRAG CALCULATIONS ////////////////
245 private double calculateFrictionDrag(Configuration configuration, FlightConditions conditions,
246 Map<RocketComponent, AerodynamicForces> map, WarningSet set) {
247 double c1 = 1.0, c2 = 1.0;
249 double mach = conditions.getMach();
254 buildCalcMap(configuration);
256 Re = conditions.getVelocity() * configuration.getLength() /
257 conditions.getAtmosphericConditions().getKinematicViscosity();
259 // System.out.printf("Re=%.3e ", Re);
261 // Calculate the skin friction coefficient (assume non-roughness limited)
262 if (configuration.getRocket().isPerfectFinish()) {
264 // System.out.printf("Perfect finish: Re=%f ",Re);
265 // Assume partial laminar layer. Roughness-limitation is checked later.
269 // System.out.printf("constant Cf=%f ",Cf);
270 } else if (Re < 5.39e5) {
272 Cf = 1.328 / MathUtil.safeSqrt(Re);
273 // System.out.printf("basic Cf=%f ",Cf);
276 Cf = 1.0 / pow2(1.50 * Math.log(Re) - 5.6) - 1700 / Re;
277 // System.out.printf("transitional Cf=%f ",Cf);
280 // Compressibility correction
283 // Below Re=1e6 no correction
286 c1 = 1 - 0.1 * pow2(mach) * (Re - 1e6) / 2e6; // transition to turbulent
288 c1 = 1 - 0.1 * pow2(mach);
295 c2 = 1 + (1.0 / Math.pow(1 + 0.045 * pow2(mach), 0.25) - 1) * (Re - 1e6) / 2e6;
297 c2 = 1.0 / Math.pow(1 + 0.045 * pow2(mach), 0.25);
302 // System.out.printf("c1=%f c2=%f\n", c1,c2);
303 // Applying continuously around Mach 1
306 } else if (mach < 1.1) {
307 Cf *= (c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2);
312 // System.out.printf("M=%f Cf=%f (smooth)\n",mach,Cf);
316 // Assume fully turbulent. Roughness-limitation is checked later.
320 // System.out.printf("LOW-TURB ");
323 Cf = 1.0 / pow2(1.50 * Math.log(Re) - 5.6);
324 // System.out.printf("NORMAL-TURB ");
327 // Compressibility correction
330 c1 = 1 - 0.1 * pow2(mach);
333 c2 = 1 / Math.pow(1 + 0.15 * pow2(mach), 0.58);
335 // Applying continuously around Mach 1
338 } else if (mach < 1.1) {
339 Cf *= c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2;
344 // System.out.printf("M=%f, Cd=%f (turbulent)\n", mach,Cf);
348 // Roughness-limited value correction term
349 double roughnessCorrection;
351 roughnessCorrection = 1 - 0.1 * pow2(mach);
352 } else if (mach > 1.1) {
353 roughnessCorrection = 1 / (1 + 0.18 * pow2(mach));
355 c1 = 1 - 0.1 * pow2(0.9);
356 c2 = 1.0 / (1 + 0.18 * pow2(1.1));
357 roughnessCorrection = c2 * (mach - 0.9) / 0.2 + c1 * (1.1 - mach) / 0.2;
360 // System.out.printf("Cf=%.3f ", Cf);
364 * Calculate the friction drag coefficient.
366 * The body wetted area is summed up and finally corrected with the rocket
367 * fineness ratio (calculated in the same iteration). The fins are corrected
368 * for thickness as we go on.
371 double finFriction = 0;
372 double bodyFriction = 0;
373 double maxR = 0, len = 0;
375 double[] roughnessLimited = new double[Finish.values().length];
376 Arrays.fill(roughnessLimited, Double.NaN);
378 for (RocketComponent c : configuration) {
380 // Consider only SymmetricComponents and FinSets:
381 if (!(c instanceof SymmetricComponent) &&
382 !(c instanceof FinSet))
385 // Calculate the roughness-limited friction coefficient
386 Finish finish = ((ExternalComponent) c).getFinish();
387 if (Double.isNaN(roughnessLimited[finish.ordinal()])) {
388 roughnessLimited[finish.ordinal()] =
389 0.032 * Math.pow(finish.getRoughnessSize() / configuration.getLength(), 0.2) *
392 // System.out.printf("roughness["+finish+"]=%.3f ",
393 // roughnessLimited[finish.ordinal()]);
397 * Actual Cf is maximum of Cf and the roughness-limited value.
398 * For perfect finish require additionally that Re > 1e6
401 if (configuration.getRocket().isPerfectFinish()) {
403 // For perfect finish require Re > 1e6
404 if ((Re > 1.0e6) && (roughnessLimited[finish.ordinal()] > Cf)) {
405 componentCf = roughnessLimited[finish.ordinal()];
406 // System.out.printf(" rl=%f Cf=%f (perfect=%b)\n",
407 // roughnessLimited[finish.ordinal()],
408 // Cf,rocket.isPerfectFinish());
410 // System.out.printf("LIMITED ");
413 // System.out.printf("NORMAL ");
418 // For fully turbulent use simple max
419 componentCf = Math.max(Cf, roughnessLimited[finish.ordinal()]);
423 // System.out.printf("compCf=%.3f ", componentCf);
428 // Calculate the friction drag:
429 if (c instanceof SymmetricComponent) {
431 SymmetricComponent s = (SymmetricComponent) c;
433 bodyFriction += componentCf * s.getComponentWetArea();
437 map.get(c).setFrictionCD(componentCf * s.getComponentWetArea()
438 / conditions.getRefArea());
441 double r = Math.max(s.getForeRadius(), s.getAftRadius());
444 len += c.getLength();
446 } else if (c instanceof FinSet) {
448 FinSet f = (FinSet) c;
449 double mac = ((FinSetCalc) calcMap.get(c)).getMACLength();
450 double cd = componentCf * (1 + 2 * f.getThickness() / mac) *
451 2 * f.getFinCount() * f.getFinArea();
455 map.get(c).setFrictionCD(cd / conditions.getRefArea());
461 // fB may be POSITIVE_INFINITY, but that's ok for us
462 double fB = (len + 0.0001) / maxR;
463 double correction = (1 + 1.0 / (2 * fB));
465 // Correct body data in map
467 for (RocketComponent c : map.keySet()) {
468 if (c instanceof SymmetricComponent) {
469 map.get(c).setFrictionCD(map.get(c).getFrictionCD() * correction);
474 // System.out.printf("\n");
475 return (finFriction + correction * bodyFriction) / conditions.getRefArea();
480 private double calculatePressureDrag(Configuration configuration, FlightConditions conditions,
481 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
483 double stagnation, base, total;
487 buildCalcMap(configuration);
489 stagnation = calculateStagnationCD(conditions.getMach());
490 base = calculateBaseCD(conditions.getMach());
493 for (RocketComponent c : configuration) {
494 if (!c.isAerodynamic())
497 // Pressure fore drag
498 double cd = calcMap.get(c).calculatePressureDragForce(conditions, stagnation, base,
503 map.get(c).setPressureCD(cd);
508 if (c instanceof SymmetricComponent) {
509 SymmetricComponent s = (SymmetricComponent) c;
511 if (radius < s.getForeRadius()) {
512 double area = Math.PI * (pow2(s.getForeRadius()) - pow2(radius));
513 cd = stagnation * area / conditions.getRefArea();
516 map.get(c).setPressureCD(map.get(c).getPressureCD() + cd);
520 radius = s.getAftRadius();
528 private double calculateBaseDrag(Configuration configuration, FlightConditions conditions,
529 Map<RocketComponent, AerodynamicForces> map, WarningSet warnings) {
533 RocketComponent prevComponent = null;
536 buildCalcMap(configuration);
538 base = calculateBaseCD(conditions.getMach());
541 for (RocketComponent c : configuration) {
542 if (!(c instanceof SymmetricComponent))
545 SymmetricComponent s = (SymmetricComponent) c;
547 if (radius > s.getForeRadius()) {
548 double area = Math.PI * (pow2(radius) - pow2(s.getForeRadius()));
549 double cd = base * area / conditions.getRefArea();
552 map.get(prevComponent).setBaseCD(cd);
556 radius = s.getAftRadius();
561 double area = Math.PI * pow2(radius);
562 double cd = base * area / conditions.getRefArea();
565 map.get(prevComponent).setBaseCD(cd);
574 public static double calculateStagnationCD(double m) {
577 pressure = 1 + pow2(m) / 4 + pow2(pow2(m)) / 40;
579 pressure = 1.84 - 0.76 / pow2(m) + 0.166 / pow2(pow2(m)) + 0.035 / pow2(m * m * m);
581 return 0.85 * pressure;
585 public static double calculateBaseCD(double m) {
587 return 0.12 + 0.13 * m * m;
595 private static final double[] axialDragPoly1, axialDragPoly2;
597 PolyInterpolator interpolator;
598 interpolator = new PolyInterpolator(
599 new double[] { 0, 17 * Math.PI / 180 },
600 new double[] { 0, 17 * Math.PI / 180 }
602 axialDragPoly1 = interpolator.interpolator(1, 1.3, 0, 0);
604 interpolator = new PolyInterpolator(
605 new double[] { 17 * Math.PI / 180, Math.PI / 2 },
606 new double[] { 17 * Math.PI / 180, Math.PI / 2 },
607 new double[] { Math.PI / 2 }
609 axialDragPoly2 = interpolator.interpolator(1.3, 0, 0, 0, 0);
614 * Calculate the axial drag from the total drag coefficient.
620 private double calculateAxialDrag(FlightConditions conditions, double cd) {
621 double aoa = MathUtil.clamp(conditions.getAOA(), 0, Math.PI);
624 // double sinaoa = conditions.getSinAOA();
625 // return cd * (1 + Math.min(sinaoa, 0.25));
628 if (aoa > Math.PI / 2)
630 if (aoa < 17 * Math.PI / 180)
631 mul = PolyInterpolator.eval(aoa, axialDragPoly1);
633 mul = PolyInterpolator.eval(aoa, axialDragPoly2);
635 if (conditions.getAOA() < Math.PI / 2)
642 private void calculateDampingMoments(Configuration configuration, FlightConditions conditions,
643 AerodynamicForces total) {
645 // Calculate pitch and yaw damping moments
646 double mul = getDampingMultiplier(configuration, conditions,
647 conditions.getPitchCenter().x);
648 double pitch = conditions.getPitchRate();
649 double yaw = conditions.getYawRate();
650 double vel = conditions.getVelocity();
652 vel = MathUtil.max(vel, 1);
654 // double Cm = total.Cm - total.CN * total.cg.x / conditions.getRefLength();
655 // System.out.printf("Damping pitch/yaw, mul=%.4f pitch rate=%.4f "+
656 // "Cm=%.4f / %.4f effect=%.4f aoa=%.4f\n", mul, pitch, total.Cm, Cm,
657 // -(mul * MathUtil.sign(pitch) * pow2(pitch/vel)),
658 // conditions.getAOA()*180/Math.PI);
660 mul *= 3; // TODO: Higher damping yields much more realistic apogee turn
662 // total.Cm -= mul * pitch / pow2(vel);
663 // total.Cyaw -= mul * yaw / pow2(vel);
664 total.setPitchDampingMoment(mul * MathUtil.sign(pitch) * pow2(pitch / vel));
665 total.setYawDampingMoment(mul * MathUtil.sign(yaw) * pow2(yaw / vel));
668 // TODO: MEDIUM: Are the rotation etc. being added correctly? sin/cos theta?
671 private double getDampingMultiplier(Configuration configuration, FlightConditions conditions,
673 if (cacheDiameter < 0) {
678 for (RocketComponent c : configuration) {
679 if (c instanceof SymmetricComponent) {
680 SymmetricComponent s = (SymmetricComponent) c;
681 area += s.getComponentPlanformArea();
682 cacheLength += s.getLength();
686 cacheDiameter = area / cacheLength;
692 mul = 0.275 * cacheDiameter / (conditions.getRefArea() * conditions.getRefLength());
693 mul *= (MathUtil.pow4(cgx) + MathUtil.pow4(cacheLength - cgx));
696 // TODO: LOW: This could be optimized a lot...
697 for (RocketComponent c : configuration) {
698 if (c instanceof FinSet) {
699 FinSet f = (FinSet) c;
700 mul += 0.6 * Math.min(f.getFinCount(), 4) * f.getFinArea() *
701 MathUtil.pow3(Math.abs(f.toAbsolute(new Coordinate(
702 ((FinSetCalc) calcMap.get(f)).getMidchordPos()))[0].x
704 (conditions.getRefArea() * conditions.getRefLength());
713 //////// The calculator map
716 protected void voidAerodynamicCache() {
717 super.voidAerodynamicCache();
725 private void buildCalcMap(Configuration configuration) {
726 Iterator<RocketComponent> iterator;
728 calcMap = new HashMap<RocketComponent, RocketComponentCalc>();
730 iterator = configuration.getRocket().iterator();
731 while (iterator.hasNext()) {
732 RocketComponent c = iterator.next();
734 if (!c.isAerodynamic())
737 calcMap.put(c, (RocketComponentCalc) Reflection.construct(BARROWMAN_PACKAGE,
738 c, BARROWMAN_SUFFIX, c));
744 public int getModID() {
745 // Only cached data is stored, return constant mod ID