diff --git a/core/src/net/sf/openrocket/aerodynamics/barrowman/RailButtonCalc.java b/core/src/net/sf/openrocket/aerodynamics/barrowman/RailButtonCalc.java index a94bd76e5..26a279fbf 100644 --- a/core/src/net/sf/openrocket/aerodynamics/barrowman/RailButtonCalc.java +++ b/core/src/net/sf/openrocket/aerodynamics/barrowman/RailButtonCalc.java @@ -58,44 +58,56 @@ public class RailButtonCalc extends RocketComponentCalc { final double notchArea = (button.getOuterDiameter() - button.getInnerDiameter()) * button.getInnerHeight(); final double refArea = outerArea - notchArea; - // accumulate Cd contribution from each rail button + // accumulate Cd contribution from each rail button. If velocity is 0 just set CDmul to a value previously + // competed for velocity MathUtil.EPSILON and skip the loop to avoid division by 0 double CDmul = 0.0; - for (int i = 0; i < button.getInstanceCount(); i++) { - - // compute boundary layer height at button location. I can't find a good reference for the - // formula, e.g. https://aerospaceengineeringblog.com/boundary-layers/ simply says it's the - // "scientific consensus". - double x = (button.toAbsolute(instanceOffsets[i]))[0].x; // location of button - double rex = calculateReynoldsNumber(x, conditions); // Reynolds number of button location - double del = 0.37 * x / Math.pow(rex, 0.2); // Boundary layer thickness - - // compute mean airspeed over button - // this assumes airspeed changes linearly through boundary layer - // and that all parts of the railbutton contribute equally to Cd, - // neither of which is true but both are plenty close enough for our purposes - - double mach; - if (buttonHt > del) { - // Case 1: button extends beyond boundary layer - // Mean velocity is 1/2 rocket velocity up to limit of boundary layer, - // full velocity after that - mach = (buttonHt - 0.5*del) * conditions.getMach()/buttonHt; - } else { - // Case 2: button is entirely within boundary layer - mach = MathUtil.map(buttonHt/2.0, 0, del, 0, conditions.getMach()); + if (conditions.getMach() > MathUtil.EPSILON) { + for (int i = 0; i < button.getInstanceCount(); i++) { + + // compute boundary layer height at button location. I can't find a good reference for the + // formula, e.g. https://aerospaceengineeringblog.com/boundary-layers/ simply says it's the + // "scientific consensus". + double x = (button.toAbsolute(instanceOffsets[i]))[0].x; // location of button + double rex = calculateReynoldsNumber(x, conditions); // Reynolds number of button location + double del = 0.37 * x / Math.pow(rex, 0.2); // Boundary layer thickness + + // compute mean airspeed over button + // this assumes airspeed changes linearly through boundary layer + // and that all parts of the railbutton contribute equally to Cd, + // neither of which is true but both are plenty close enough for our purposes + + double mach; + if (buttonHt > del) { + // Case 1: button extends beyond boundary layer + // Mean velocity is 1/2 rocket velocity up to limit of boundary layer, + // full velocity after that + mach = (buttonHt - 0.5*del) * conditions.getMach()/buttonHt; + } else { + // Case 2: button is entirely within boundary layer + mach = MathUtil.map(buttonHt/2.0, 0, del, 0, conditions.getMach()); + } + + // look up Cd as function of speed. It's pretty constant as a function of Reynolds + // number when slow, so we can just use a function of Mach number + double cd = MathUtil.interpolate(cdDomain, cdRange, mach); + + // Since later drag force calculations don't consider boundary layer, compute "effective Cd" + // based on rocket velocity + cd = cd * MathUtil.pow2(mach)/MathUtil.pow2(conditions.getMach()); + + // add to CDmul + CDmul += cd; + } - // look up Cd as function of speed. It's pretty constant as a function of Reynolds - // number when slow, so we can just use a function of Mach number - double cd = MathUtil.interpolate(cdDomain, cdRange, mach); - - // Since later drag force calculations don't consider boundary layer, compute "effective Cd" - // based on rocket velocity - cd = cd * MathUtil.pow2(mach)/MathUtil.pow2(conditions.getMach()); + // since we'll be multiplying by the instance count up in BarrowmanCalculator, + // we want to return the mean CD instead of the total + CDmul /= button.getInstanceCount(); - // add to CDmul - CDmul += cd; - } + } else { + // value at velocity of MathUtil.EPSILON + CDmul = 8.786395072609939E-4; + } return CDmul * stagnationCD * refArea / conditions.getRefArea(); } diff --git a/core/test/net/sf/openrocket/aerodynamics/BarrowmanCalculatorTest.java b/core/test/net/sf/openrocket/aerodynamics/BarrowmanCalculatorTest.java index f51b9ee99..a3ec7c07b 100644 --- a/core/test/net/sf/openrocket/aerodynamics/BarrowmanCalculatorTest.java +++ b/core/test/net/sf/openrocket/aerodynamics/BarrowmanCalculatorTest.java @@ -20,11 +20,13 @@ import net.sf.openrocket.rocketcomponent.FlightConfiguration; import net.sf.openrocket.rocketcomponent.NoseCone; import net.sf.openrocket.rocketcomponent.ParallelStage; import net.sf.openrocket.rocketcomponent.PodSet; +import net.sf.openrocket.rocketcomponent.RailButton; import net.sf.openrocket.rocketcomponent.Rocket; import net.sf.openrocket.rocketcomponent.Transition; import net.sf.openrocket.rocketcomponent.TrapezoidFinSet; import net.sf.openrocket.startup.Application; import net.sf.openrocket.util.Coordinate; +import net.sf.openrocket.util.MathUtil; import net.sf.openrocket.util.TestRockets; public class BarrowmanCalculatorTest { @@ -496,5 +498,70 @@ public class BarrowmanCalculatorTest { testCP = testCalc.getCP(testConfig, testConditions, warnings).x; assertEquals("should be warning from podset airframe overlap", 1, warnings.size()); } + + /** + * Tests railbutton drag. Really is testing instancing more than actual drag calculations, and making + * sure we don't divide by 0 when not moving + */ + @Test + public void testRailButtonDrag() { + // minimal rocket with nothing on it but two railbuttons + final Rocket rocket = new Rocket(); + final AxialStage stage = new AxialStage(); + rocket.addChild(stage); + + // phantom tubes have no drag to confuse things + final BodyTube phantom = new BodyTube(); + phantom.setOuterRadius(0); + stage.addChild(phantom); + + // set up test environment + WarningSet warnings = new WarningSet(); + final FlightConfiguration config = rocket.getSelectedConfiguration(); + final FlightConditions conditions = new FlightConditions(config); + final BarrowmanCalculator calc = new BarrowmanCalculator(); + + // part 1: instancing + + // Put two individual railbuttons and get their CD + final RailButton button1 = new RailButton(); + button1.setInstanceCount(1); + button1.setAxialOffset(1.0); + phantom.addChild(button1); + + final RailButton button2 = new RailButton(); + button2.setInstanceCount(1); + button2.setAxialOffset(2.0); + phantom.addChild(button2); + + final AerodynamicForces individualForces = calc.getAerodynamicForces(config, conditions, warnings); + final double individualCD = individualForces.getCD(); + + // get rid of individual buttons and put in a railbutton set with two instances at same locations as original + // railbuttons + phantom.removeChild(button1); + phantom.removeChild(button2); + + final RailButton buttons = new RailButton(); + buttons.setInstanceCount(2); + buttons.setAxialOffset(1.0); + buttons.setInstanceSeparation(1.0); + + final AerodynamicForces pairForces = calc.getAerodynamicForces(config, conditions, warnings); + final double pairCD = pairForces.getCD(); + + assertEquals("two individual railbuttons should have same CD as a pair", individualCD, pairCD, EPSILON); + + // part 2: test at Mach 0 + conditions.setMach(MathUtil.EPSILON); + final AerodynamicForces epsForces = calc.getAerodynamicForces(config, conditions, warnings); + final double epsCD = epsForces.getCD(); + + conditions.setMach(0); + final AerodynamicForces zeroForces = calc.getAerodynamicForces(config, conditions, warnings); + final double zeroCD = zeroForces.getCD(); + assertEquals("drag at mach 0 should equal drag at mach MathUtil.EPSILON", epsCD, zeroCD, EPSILON); + } } +