Merge pull request #2066 from JoePfeiffer/fix-tubefins

Fix interstices area calculation and update pressure drag calculation for tube fins

core/src/net/sf/openrocket/aerodynamics/barrowman/TubeCalc.java

@@ -14,73 +14,64 @@ public abstract class TubeCalc extends RocketComponentCalc {
 	
 	private final static Logger log = LoggerFactory.getLogger(TubeFinSetCalc.class);
 
-	// air density (standard conditions)
+	private final Tube tube;
-	private final double rho = 1.225; // kg/m^3
-	
 	private final double diameter;
 	private final double length;
 	protected final double innerArea;
 	private final double totalArea;
 	private final double frontalArea;
+	private final double epsilon;
 	
 	public TubeCalc(RocketComponent component) {
 		super(component);
 
-		Tube tube = (Tube)component;
+		tube = (Tube)component;
 		
 		length = tube.getLength();
 		diameter = 2 * tube.getInnerRadius();
 		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
+			// Current atmospheric conditions
-			final double T = conditions.getAtmosphericConditions().getTemperature();
+			final double p = conditions.getAtmosphericConditions().getPressure();
-			
+			final double t = conditions.getAtmosphericConditions().getTemperature();
-			// Volume flow rate (t)
+			final double rho = conditions.getAtmosphericConditions().getDensity();
-			final double Q = conditions.getVelocity() * innerArea;
+			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 = conditions.getVelocity() * diameter /
+			final double Re = v * diameter / conditions.getAtmosphericConditions().getKinematicViscosity();
-				conditions.getAtmosphericConditions().getKinematicViscosity();
 			
-			// friction coefficient (for smooth tube interior) (e)
+			// friction coefficient using Swamee-Jain equation
-			final double lambda = 1/MathUtil.pow2(2 * Math.log(0.5625 * Math.pow(Re, 0.875)) - 0.8);
+			double f = 0.25/MathUtil.pow2(Math.log10((epsilon / (3.7 * diameter) + 5.74/Math.pow(Re, 0.9))));
 
-			// pressure drop (e)
+			// If we're supersonic, apply a correction
-			final double P0 = 100; // standard pressure
+			if (conditions.getMach() > 1) {
-			final double T0 = 273.15; // standard temperature
+			 	f = f / conditions.getBeta();
-			deltap = (lambda * 8 * length * rho * MathUtil.pow2(Q) * T * P0) /
+			}
-				(MathUtil.pow2(Math.PI) * Math.pow(diameter, 5) * T0 * conditions.getAtmosphericConditions().getPressure());
+
-		} else {
+			// pressure drop using Darcy-Weissbach equation
-			deltap = 0.0;
+			deltap = f * (length * rho * MathUtil.pow2(v)) / (2 * diameter);
+
+			// drag coefficient of tube interior from pressure drop
+			tubeCD = 2 * (deltap * innerArea) / (rho * MathUtil.pow2(v) * innerArea);
 		}
 		   
 		// convert to CD and return
-		return (deltap * innerArea + 0.7 * stagnationCD * frontalArea) / conditions.getRefArea();
+		final double cd = (tubeCD * innerArea + 0.7*(stagnationCD + baseCD) * frontalArea) / conditions.getRefArea();
+		return cd;
 	}
 }

core/src/net/sf/openrocket/aerodynamics/barrowman/TubeFinSetCalc.java

@@ -26,6 +26,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
 	private final double bodyRadius;
@@ -51,7 +53,7 @@ public class TubeFinSetCalc extends TubeCalc {
 			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);
@@ -76,36 +78,52 @@ public class TubeFinSetCalc extends TubeCalc {
 		// aspect ratio.
 		ar = 2 * innerRadius / chord;
 
+
+		// Some trigonometry...
+		// We need a triangle with the following three sides:
+		// d is from the center of the body tube to a tangent point on the tube fin
+		// outerRadius is from the center of the tube fin to the tangent point.  Note that
+		//     d and outerRadius are at right angles
+		// bodyRadius + outerRadius is from the center of the body tube to the center of the tube fin.
+		//     This is the hypotenuse of the right triangle.
+		
+		// Find length of d
+		final double d = Math.sqrt(MathUtil.pow2(bodyRadius + outerRadius) - MathUtil.pow2(outerRadius));
+
+		// 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 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;
+		
+		// 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;
+		
 		// wetted area for friction drag calculation.  We don't consider the inner surface of the tube;
 		// that affects the pressure drop through the tube and so (indirecctly) affects the pressure drag.
 			
-		// Area of the outer surface of tubes.  Since roughly half
+		// Area of the outer surface of a tube, not including portion masked by interstice
-		// of the area is "masked" by the interstices between the tubes and the
+		final double outerArea = chord * 2.0 * (Math.PI - theta1) * outerRadius;
-		// body tube, only consider the other half of the area (so only multiplying by pi instead of 2*pi)
-		final double outerArea = chord * Math.PI * outerRadius;
 			
-		// Surface area of the portion of the body tube masked by the tube fins, per tube
+		// Surface area of the portion of the body tube masked by the tube fin.  We'll subtract it from
-		final BodyTube parent = (BodyTube) tubes.getParent();
+		// the tube fin area rather than go in and change the body tube surface area calculation. If tube
-		final double maskedArea = chord * 2.0 * Math.PI * bodyRadius / tubeCount;
+		// 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;
-		log.debug("wetted area of tube fins " + wettedArea);
-
-		// frontal area of interstices between tubes for pressure drag calculation.
-		// We'll treat them as a closed blunt object.
-
-		// area of disk passing through tube fin centers
-		final double tubeDiskArea = Math.PI * MathUtil.pow2(bodyRadius + outerRadius);
-
-		// half of combined area of tube fin exteriors.  Deliberately using the outer radius here since we
-		// calculate pressure drag from the tube walls in TubeCalc
-		final double tubeOuterArea = tubeCount * Math.PI * MathUtil.pow2(outerRadius) / 2.0;
-
-		// body tube area
-		final double bodyTubeArea = Math.PI * MathUtil.pow2(bodyRadius);
-
-		// area of an interstice
-		intersticeArea = (tubeDiskArea - tubeOuterArea - bodyTubeArea) / tubeCount;
 
 		// Precompute most of CNa.  Equation comes from Ribner, "The ring airfoil in nonaxial
 		// flow", Journal of the Aeronautical Sciences 14(9) pp 529-530 (1947) equation (5).
@@ -246,8 +264,6 @@ public class TubeFinSetCalc extends TubeCalc {
 
 	@Override
 	public double calculateFrictionCD(FlightConditions conditions, double componentCf, WarningSet warnings) {
-		warnings.addAll(geometryWarnings);
-
 		final double frictionCD = componentCf * wettedArea / conditions.getRefArea();
 
 		return frictionCD;
@@ -258,18 +274,10 @@ public class TubeFinSetCalc extends TubeCalc {
 					  double stagnationCD, double baseCD, WarningSet warnings) {
 		
 	    warnings.addAll(geometryWarnings);
+						    
 		final double cd = super.calculatePressureCD(conditions, stagnationCD, baseCD, warnings) +
 					(stagnationCD + baseCD) * intersticeArea / conditions.getRefArea();
 
 	    return cd;
 	}
-	
-	private static int calculateInterferenceFinCount(TubeFinSet component) {
-		RocketComponent parent = component.getParent();
-		if (parent == null) {
-			throw new IllegalStateException("fin set without parent component");
-		}
-		
-		return 3 * component.getFinCount();
-	}
 }

core/src/net/sf/openrocket/rocketcomponent/BodyTube.java

@@ -550,4 +550,19 @@ public class BodyTube extends SymmetricComponent implements BoxBounded, MotorMou
 		// The motor config also has listeners, so clear them as well
 		getDefaultMotorConfig().clearConfigListeners();
 	}
+
+	/**
+	 * The first time we add a TubeFinSet to the component tree, inherit the tube thickness from
+	 * the parent body tube
+	 */
+	@Override
+	public final void addChild(RocketComponent component, int index, boolean trackStage) {
+		super.addChild(component, index, trackStage);
+		if (component instanceof TubeFinSet) {
+			TubeFinSet finset = (TubeFinSet) component;
+			if (Double.isNaN(finset.getThickness())) {
+				finset.setThickness(getThickness());
+			}
+		}
+	}
 }

core/src/net/sf/openrocket/rocketcomponent/TubeFinSet.java

@@ -24,7 +24,7 @@ public class TubeFinSet extends Tube implements AxialPositionable, BoxBounded, R
 	
 	private boolean autoRadius = true; // Radius chosen automatically based on parent component
 	private double outerRadius = DEFAULT_RADIUS;
-	protected double thickness = 0.002;
+	protected double thickness = Double.NaN;
 	private AngleMethod angleMethod = AngleMethod.FIXED;
 	protected RadiusMethod radiusMethod = RadiusMethod.RELATIVE;
 
@@ -49,7 +49,7 @@ public class TubeFinSet extends Tube implements AxialPositionable, BoxBounded, R
 	
 	
 	/**
-	 * New FinSet with given number of fins and given base rotation angle.
+	 * New TubeFinSet with default values
 	 * Sets the component relative position to POSITION_RELATIVE_BOTTOM,
 	 * i.e. fins are positioned at the bottom of the parent component.
 	 */
@@ -146,6 +146,7 @@ public class TubeFinSet extends Tube implements AxialPositionable, BoxBounded, R
 	 * Sets whether the radius is selected automatically or not.
 	 */
 	public void setOuterRadiusAutomatic(boolean auto) {
+		
 		for (RocketComponent listener : configListeners) {
 			if (listener instanceof TubeFinSet) {
 				((TubeFinSet) listener).setOuterRadiusAutomatic(auto);
@@ -195,8 +196,9 @@ public class TubeFinSet extends Tube implements AxialPositionable, BoxBounded, R
 
 		if ((this.thickness == thickness))
 			return;
+
 		this.thickness = MathUtil.clamp(thickness, 0, getOuterRadius());
-		fireComponentChangeEvent(ComponentChangeEvent.MASS_CHANGE);
+		fireComponentChangeEvent(ComponentChangeEvent.BOTH_CHANGE);
 		clearPreset();
 	}