diff --git a/core/src/net/sf/openrocket/aerodynamics/barrowman/TubeCalc.java b/core/src/net/sf/openrocket/aerodynamics/barrowman/TubeCalc.java index 82e6e4662..ed34e7476 100644 --- a/core/src/net/sf/openrocket/aerodynamics/barrowman/TubeCalc.java +++ b/core/src/net/sf/openrocket/aerodynamics/barrowman/TubeCalc.java @@ -14,15 +14,13 @@ public abstract class TubeCalc extends RocketComponentCalc { private final static Logger log = LoggerFactory.getLogger(TubeFinSetCalc.class); - // air density (standard conditions) - private final double rho = 1.225; // kg/m^3 - + private final Tube tube; private final double diameter; private final double length; protected final double innerArea; private final double totalArea; private final double frontalArea; - private final Tube tube; + private final double epsilon; public TubeCalc(RocketComponent component) { super(component); @@ -34,63 +32,49 @@ public abstract class TubeCalc extends RocketComponentCalc { innerArea = Math.PI * MathUtil.pow2(tube.getInnerRadius()); totalArea = Math.PI * MathUtil.pow2(tube.getOuterRadius()); frontalArea = totalArea - innerArea; + epsilon = tube.getFinish().getRoughnessSize(); // roughness; note we don't maintain surface roughness of + // interior separately from exterior. } @Override public double calculatePressureCD(FlightConditions conditions, double stagnationCD, double baseCD, WarningSet warnings) { - - // These calculations come from a mix of theoretical and empirical - // results, and are marked with (t) for theoretical and (e) for empirical. - // The theoretical results should not be modified; the empirical can be adjusted - // to better simulate real rockets as we get data. - - // For the sources of the empirical formulas, see Carello, Ivanov, and Mazza, - // "Pressure drop in pipe lines for compressed air: comparison between experimental - // and theoretical analysis", Transactions on Engineering Sciences vol 18, - // ISSN 1743-35331998, 1998. - - // For the rockets for which we have data, the effect of the stagnation CD appears to be - // overstated. This code multiplies it be a factor of 0.7 to better match experimental - // data // Need to check for tube inner area 0 in case of rockets using launch lugs with - // an inner radius of 0 to emulate rail buttons (or just weird rockets, of course) - + // an inner radius of 0 to emulate rail guides (or just weird rockets, of course) + double tubeCD = 0.0; double deltap; if (innerArea > MathUtil.EPSILON) { - // Temperature and Pressure - final double T = conditions.getAtmosphericConditions().getTemperature(); - final double P = conditions.getAtmosphericConditions().getPressure(); - - // Volume flow rate (t) - final double Q = conditions.getVelocity() * innerArea; - - // Air viscosity - final double mu = conditions.getAtmosphericConditions().getKinematicViscosity(); - - // Air density - final double rho = 1.225; // at standard temperature and pressure + // Current atmospheric conditions + final double p = conditions.getAtmosphericConditions().getPressure(); + final double t = conditions.getAtmosphericConditions().getTemperature(); + final double rho = conditions.getAtmosphericConditions().getDensity(); + final double v = conditions.getVelocity(); // Reynolds number (note Reynolds number for the interior of a pipe is based on diameter, // not length (t)) - final double Re = (4.0 * rho * Q) / (Math.PI * diameter * mu); + final double Re = v * diameter / conditions.getAtmosphericConditions().getKinematicViscosity(); - // friction coefficient (for smooth tube interior) (e) - final double lambda = 1/MathUtil.pow2(2 * Math.log(0.5625 * Math.pow(Re, 0.875)) - 0.8); - - // pressure drop (e) - final double P0 = 101325; // standard pressure - final double T0 = 273.15; // standard temperature - deltap = ((lambda * 8 * length * rho * MathUtil.pow2(Q)) / (MathUtil.pow2(Math.PI) * Math.pow(diameter, 5)) * (T/T0) * (P0/P)); - } else { - deltap = 0.0; + // friction coefficient using Swamee-Jain equation + double f = 0.25/MathUtil.pow2(Math.log10((epsilon / (3.7 * diameter) + 5.74/Math.pow(Re, 0.9)))); + + // If we're supersonic, apply a correction + // if (conditions.getMach() > 1) { + // f = f / conditions.getBeta(); + // } + + // pressure drop using Darcy-Weissbach equation + deltap = f * (length * rho * MathUtil.pow2(v)) / (2 * diameter); + System.out.println(tube + ", v " + v + ", Re " + Re + ", p " + p + ": " + "deltap " + deltap); + + // drag coefficient of tube interior from pressure drop + tubeCD = 2 * (deltap * innerArea) / (rho * MathUtil.pow2(v) * innerArea); } // convert to CD and return - final double cdpress = 2.0 * deltap / (conditions.getAtmosphericConditions().getDensity() * MathUtil.pow2(conditions.getVelocity())); - final double cd = (cdpress * innerArea + 0.43*(stagnationCD + baseCD) * frontalArea)/conditions.getRefArea(); - + System.out.println(tube + " tube CD " + tubeCD + ", stagnationCD " + stagnationCD + ", baseCD " + baseCD + ", inner area " + innerArea + ", frontal area " + frontalArea); + final double cd = (tubeCD * innerArea + (stagnationCD + baseCD) * frontalArea) / conditions.getRefArea(); + System.out.println(tube + " cd " + cd); return cd; } } diff --git a/core/src/net/sf/openrocket/aerodynamics/barrowman/TubeFinSetCalc.java b/core/src/net/sf/openrocket/aerodynamics/barrowman/TubeFinSetCalc.java index d0baed225..1a2e8f97e 100644 --- a/core/src/net/sf/openrocket/aerodynamics/barrowman/TubeFinSetCalc.java +++ b/core/src/net/sf/openrocket/aerodynamics/barrowman/TubeFinSetCalc.java @@ -25,6 +25,8 @@ public class TubeFinSetCalc extends TubeCalc { private static final double STALL_ANGLE = (20 * Math.PI / 180); private final double[] poly = new double[6]; + + private final TubeFinSet tubes; // parameters straight from configuration; we'll be grabbing them once // so code is a bit shorter elsewhere @@ -44,14 +46,14 @@ public class TubeFinSetCalc extends TubeCalc { private final double cnaconst; protected final WarningSet geometryWarnings = new WarningSet(); - + public TubeFinSetCalc(RocketComponent component) { super(component); if (!(component instanceof TubeFinSet)) { throw new IllegalArgumentException("Illegal component type " + component); } - final TubeFinSet tubes = (TubeFinSet) component; + tubes = (TubeFinSet) component; if (tubes.getTubeSeparation() > MathUtil.EPSILON) { geometryWarnings.add(Warning.TUBE_SEPARATION); @@ -88,24 +90,26 @@ public class TubeFinSetCalc extends TubeCalc { // Find length of d final double d = Math.sqrt(MathUtil.pow2(bodyRadius + outerRadius) - MathUtil.pow2(outerRadius)); - // Area of diamond consisting of triangle reflected on its hypotenuse + // Area of diamond formed by mirroring triangle on its hypotenuse (same area as rectangle + // formed by d and outerarea, but it *isn't* that rectangle) double a = d * outerRadius; - + // angle between outerRadius and bodyRadius+outerRadius final double theta1 = Math.acos(outerRadius/(outerRadius + bodyRadius)); - // area of arc from tube fin, doubled so we have area to remove from diamond + // area of arc from tube fin, doubled to get both halves of diamond final double a1 = MathUtil.pow2(outerRadius) * theta1; - + // angle between bodyRadius+outerRadius and d final double theta2 = Math.PI/2.0 - theta1; + System.out.println("theta2 " + theta2); // area of arc from body tube. Doubled so we have area to remove from diamond final double a2 = MathUtil.pow2(bodyRadius) * theta2; - + // area of interstice for one tube fin intersticeArea = (a - a1 - a2); - + // for comparison, what's the area of a tube fin? double tubeArea = MathUtil.pow2(outerRadius) * Math.PI; @@ -113,13 +117,18 @@ public class TubeFinSetCalc extends TubeCalc { // that affects the pressure drop through the tube and so (indirecctly) affects the pressure drag. // Area of the outer surface of a tube, not including portion masked by interstice - final double outerArea = chord * 2 * (Math.PI - theta1) * outerRadius; + final double outerArea = chord * 2.0 * (Math.PI - theta1) * outerRadius; + + // Area of inner surface of a tube + final double innerArea = chord * 2.0 * Math.PI * innerRadius; - // Surface area of the portion of the body tube masked by the tube fins, per tube - final BodyTube parent = (BodyTube) tubes.getParent(); - final double maskedArea = chord * 2.0 * Math.PI * bodyRadius / tubeCount; + // Surface area of the portion of the body tube masked by the tube fin. We'll subtract it from + // the tube fin area rather than go in and change the body tube surface area calculation. If tube + // fin and body tube roughness aren't the same this will result in an inaccuracy. + final double maskedArea = chord * 2.0 * theta2 * bodyRadius; - wettedArea = outerArea - maskedArea; + wettedArea = innerArea + outerArea - maskedArea; + System.out.println(tubes + " outer " + outerArea + ", masked " + maskedArea); log.debug("wetted area of tube fin " + wettedArea); // Precompute most of CNa. Equation comes from Ribner, "The ring airfoil in nonaxial @@ -261,6 +270,7 @@ public class TubeFinSetCalc extends TubeCalc { @Override public double calculateFrictionCD(FlightConditions conditions, double componentCf, WarningSet warnings) { + System.out.println(tubes + "wetted area " + wettedArea); final double frictionCD = componentCf * wettedArea / conditions.getRefArea(); return frictionCD; @@ -271,9 +281,10 @@ public class TubeFinSetCalc extends TubeCalc { double stagnationCD, double baseCD, WarningSet warnings) { warnings.addAll(geometryWarnings); + System.out.println(tubes + " stag CD " + stagnationCD + ", base CD " + baseCD); final double cd = super.calculatePressureCD(conditions, stagnationCD, baseCD, warnings) + - (stagnationCD + baseCD) * intersticeArea / conditions.getRefArea(); + (stagnationCD + baseCD) * intersticeArea / conditions.getRefArea(); return cd; }