Create Tube parent class for LaunchLug and TubeFinSet, for use with TubeCalc
Create TubeCalc parent class for LaunchLugCalc and TubeFinSetCalc Replace tube drag calculation with new code, no longer using the three-fin approximation from 15.03
This commit is contained in:
JoePfeiffer
parent
57d6d44c4f
commit
265a6e0b7b
@ -585,6 +585,7 @@ public class BarrowmanCalculator extends AbstractAerodynamicCalculator {
|
||||
final InstanceMap imap = configuration.getActiveInstances();
|
||||
for(Map.Entry<RocketComponent, ArrayList<InstanceContext>> entry: imap.entrySet() ) {
|
||||
final RocketComponent c = entry.getKey();
|
||||
log.debug("component " + c);
|
||||
|
||||
if (!c.isAerodynamic())
|
||||
continue;
|
||||
@ -593,7 +594,7 @@ public class BarrowmanCalculator extends AbstractAerodynamicCalculator {
|
||||
final ArrayList<InstanceContext> contextList = entry.getValue();
|
||||
for(InstanceContext context: contextList ) {
|
||||
|
||||
// Pressure fore drag
|
||||
// Pressure drag
|
||||
double cd = calcMap.get(c).calculatePressureCD(conditions, stagnation, base,
|
||||
warningSet);
|
||||
total += cd;
|
||||
|
||||
@ -10,18 +10,8 @@ import net.sf.openrocket.util.Transformation;
|
||||
|
||||
public class LaunchLugCalc extends TubeCalc {
|
||||
|
||||
private final double CDmul;
|
||||
private final double refArea;
|
||||
|
||||
public LaunchLugCalc(RocketComponent component) {
|
||||
super(component);
|
||||
|
||||
LaunchLug lug = (LaunchLug)component;
|
||||
double ld = lug.getLength() / (2*lug.getOuterRadius());
|
||||
|
||||
CDmul = Math.max(1.3 - ld, 1);
|
||||
refArea = Math.PI * MathUtil.pow2(lug.getOuterRadius()) -
|
||||
Math.PI * MathUtil.pow2(lug.getInnerRadius()) * Math.max(1 - ld, 0);
|
||||
}
|
||||
|
||||
@Override
|
||||
@ -35,12 +25,4 @@ public class LaunchLugCalc extends TubeCalc {
|
||||
// launch lug doesn't add enough area to worry about
|
||||
return 0;
|
||||
}
|
||||
|
||||
@Override
|
||||
public double calculatePressureCD(FlightConditions conditions,
|
||||
double stagnationCD, double baseCD, WarningSet warnings) {
|
||||
|
||||
return CDmul*stagnationCD * refArea / conditions.getRefArea();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -10,19 +10,38 @@ import net.sf.openrocket.util.Transformation;
|
||||
|
||||
public abstract class TubeCalc extends RocketComponentCalc {
|
||||
|
||||
private final double CDmul;
|
||||
protected final double refArea;
|
||||
private final double diameter;
|
||||
private final double length;
|
||||
protected double refArea;
|
||||
|
||||
public TubeCalc(RocketComponent component) {
|
||||
super(component);
|
||||
|
||||
Tube tube = (Tube)component;
|
||||
|
||||
double ld = tube.getLength() / (2*tube.getOuterRadius());
|
||||
|
||||
CDmul = Math.max(1.3 - ld, 1);
|
||||
refArea = Math.PI * MathUtil.pow2(tube.getOuterRadius()) -
|
||||
Math.PI * MathUtil.pow2(tube.getInnerRadius()) * Math.max(1 - ld, 0);
|
||||
length = tube.getLength();
|
||||
diameter = 2 * tube.getInnerRadius();
|
||||
refArea = Math.PI * MathUtil.pow2(tube.getInnerRadius());
|
||||
}
|
||||
|
||||
@Override
|
||||
public double calculatePressureCD(FlightConditions conditions,
|
||||
double stagnationCD, double baseCD, WarningSet warnings) {
|
||||
|
||||
// calculation of pressure drop through pipe from "Atlas Copco Air Compendium", 1975,
|
||||
// quoted as equation 14 in 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.
|
||||
|
||||
// Volume flow rate
|
||||
final double Q = conditions.getVelocity() * refArea;
|
||||
|
||||
// pressure drop
|
||||
final double deltap = 1.6 * 1000 * Math.pow(Q, 1.85) * Math.pow(Q, 1.85) * length /
|
||||
(Math.pow(diameter, 5) * conditions.getAtmosphericConditions().getPressure());
|
||||
|
||||
// convert to CD and return
|
||||
return deltap * refArea / conditions.getRefArea();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
@ -29,14 +29,14 @@ import org.slf4j.LoggerFactory;
|
||||
* Uses a complete clone of FinSetCalc modelling each tube fin as 3 individual fins. It does not correctly account for
|
||||
* fin & tube fin interference.
|
||||
*
|
||||
* Changes to BarrowmanCalculator's calculateFrictionDrag are also probably required.
|
||||
*
|
||||
* @author kruland
|
||||
*
|
||||
*/
|
||||
public class TubeFinSetCalc extends TubeCalc {
|
||||
|
||||
private final static Logger log = LoggerFactory.getLogger(FinSetCalc.class);
|
||||
private final static Logger log = LoggerFactory.getLogger(TubeFinSetCalc.class);
|
||||
|
||||
final double intersticeArea;
|
||||
|
||||
private static final double STALL_ANGLE = (20 * Math.PI / 180);
|
||||
|
||||
@ -63,7 +63,6 @@ public class TubeFinSetCalc extends TubeCalc {
|
||||
private final double[] poly = new double[6];
|
||||
|
||||
private final double wettedArea;
|
||||
private final double frontalArea;
|
||||
|
||||
private final double thickness;
|
||||
private final double bodyRadius;
|
||||
@ -89,27 +88,28 @@ public class TubeFinSetCalc extends TubeCalc {
|
||||
geometryWarnings.add(Warning.TUBE_OVERLAP);
|
||||
}
|
||||
|
||||
// precompute geometry
|
||||
// precompute geometry. This will be the geometry of a single tube, since BarrowmanCalculator
|
||||
// iterates across them. Doesn't consider interference between them; that should only be relevant for
|
||||
// fins that are either separated or overlapping.
|
||||
bodyRadius = tubes.getBodyRadius();
|
||||
|
||||
// wetted area for friction drag calculation
|
||||
// Area of inner surface of tubes.
|
||||
final double innerArea = tubes.getFinCount() * tubes.getLength() * 2.0 * Math.PI * tubes.getInnerRadius();
|
||||
// 1. 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
|
||||
// of the area is "masked" by the interstices between the tubes and the
|
||||
// body tube, only consider half the area
|
||||
final double outerArea = tubes.getFinCount() * tubes.getLength() * Math.PI * tubes.getOuterRadius();
|
||||
// body tube, only consider the other half of the area
|
||||
final double outerArea = tubes.getLength() * Math.PI * tubes.getOuterRadius();
|
||||
|
||||
// Area of the portion of the body tube masked by the tube fins
|
||||
// Surface area of the portion of the body tube masked by the tube fins
|
||||
final BodyTube parent = (BodyTube) tubes.getParent();
|
||||
final double maskedArea = tubes.getLength() * 2.0 * Math.PI * parent.getOuterRadius();
|
||||
|
||||
wettedArea = innerArea + outerArea - maskedArea;
|
||||
wettedArea = outerArea - maskedArea;
|
||||
log.debug("wetted area of tube fins " + wettedArea);
|
||||
|
||||
// effective frontal area for pressure drag calculation
|
||||
// model interstices between tube fins and body tube as closed blunt area
|
||||
// 2. 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 + tubes.getOuterRadius());
|
||||
@ -120,9 +120,9 @@ public class TubeFinSetCalc extends TubeCalc {
|
||||
// body tube area
|
||||
final double bodyTubeArea = Math.PI * MathUtil.pow2(bodyRadius);
|
||||
|
||||
// refArea is effective frontal area of a tube, from LaunchLugCalc
|
||||
frontalArea = tubes.getFinCount() * refArea + tubeDiskArea - tubeInnerArea - bodyTubeArea;
|
||||
log.debug("effective frontal area " + frontalArea);
|
||||
// area of an interstice
|
||||
intersticeArea = (tubeDiskArea - tubeInnerArea - bodyTubeArea) / tubes.getFinCount();
|
||||
log.debug("interstice area " + intersticeArea);
|
||||
|
||||
thickness = fin.getThickness();
|
||||
finCount = 3 * fin.getFinCount();
|
||||
@ -643,18 +643,18 @@ public class TubeFinSetCalc extends TubeCalc {
|
||||
log.debug("frictionCD " + frictionCD);
|
||||
return frictionCD;
|
||||
}
|
||||
|
||||
|
||||
@Override
|
||||
public double calculatePressureCD(FlightConditions conditions,
|
||||
double stagnationCD, double baseCD, WarningSet warnings) {
|
||||
|
||||
warnings.addAll(geometryWarnings);
|
||||
|
||||
final double cd = (stagnationCD + baseCD) * frontalArea / conditions.getRefArea();
|
||||
final double cd = super.calculatePressureCD(conditions, stagnationCD, baseCD, warnings) +
|
||||
(stagnationCD + baseCD) * intersticeArea / conditions.getRefArea();
|
||||
log.debug("pressure CD " + cd);
|
||||
|
||||
return cd;
|
||||
}
|
||||
}
|
||||
|
||||
private static int calculateInterferenceFinCount(TubeFinSet component) {
|
||||
RocketComponent parent = component.getParent();
|
||||
|
||||
@ -0,0 +1,710 @@
|
||||
package net.sf.openrocket.aerodynamics.barrowman;
|
||||
|
||||
import static java.lang.Math.pow;
|
||||
import static net.sf.openrocket.util.MathUtil.pow2;
|
||||
|
||||
import java.util.Arrays;
|
||||
|
||||
import net.sf.openrocket.aerodynamics.AerodynamicForces;
|
||||
import net.sf.openrocket.aerodynamics.FlightConditions;
|
||||
import net.sf.openrocket.aerodynamics.Warning;
|
||||
import net.sf.openrocket.aerodynamics.Warning.Other;
|
||||
import net.sf.openrocket.aerodynamics.WarningSet;
|
||||
import net.sf.openrocket.rocketcomponent.BodyTube;
|
||||
import net.sf.openrocket.rocketcomponent.FinSet;
|
||||
import net.sf.openrocket.rocketcomponent.RocketComponent;
|
||||
import net.sf.openrocket.rocketcomponent.TubeFinSet;
|
||||
import net.sf.openrocket.util.Coordinate;
|
||||
import net.sf.openrocket.util.LinearInterpolator;
|
||||
import net.sf.openrocket.util.MathUtil;
|
||||
import net.sf.openrocket.util.PolyInterpolator;
|
||||
import net.sf.openrocket.util.Transformation;
|
||||
|
||||
import org.slf4j.Logger;
|
||||
import org.slf4j.LoggerFactory;
|
||||
|
||||
/**
|
||||
* Preliminary computation of tube fin aerodynamics.
|
||||
*
|
||||
* Uses a complete clone of FinSetCalc modelling each tube fin as 3 individual fins. It does not correctly account for
|
||||
* fin & tube fin interference.
|
||||
*
|
||||
* Changes to BarrowmanCalculator's calculateFrictionDrag are also probably required.
|
||||
*
|
||||
* @author kruland
|
||||
*
|
||||
*/
|
||||
public class TubeFinSetCalc extends LaunchLugCalc {
|
||||
|
||||
private final static Logger log = LoggerFactory.getLogger(FinSetCalc.class);
|
||||
|
||||
private static final double STALL_ANGLE = (20 * Math.PI / 180);
|
||||
|
||||
/** Number of divisions in the fin chords. */
|
||||
protected static final int DIVISIONS = 48;
|
||||
|
||||
protected double macLength = Double.NaN; // MAC length
|
||||
protected double macLead = Double.NaN; // MAC leading edge position
|
||||
protected double macSpan = Double.NaN; // MAC spanwise position
|
||||
protected double finArea = Double.NaN; // Fin area
|
||||
protected double ar = Double.NaN; // Fin aspect ratio
|
||||
protected double span = Double.NaN; // Fin span
|
||||
protected double cosGamma = Double.NaN; // Cosine of midchord sweep angle
|
||||
protected double cosGammaLead = Double.NaN; // Cosine of leading edge sweep angle
|
||||
protected double rollSum = Double.NaN; // Roll damping sum term
|
||||
|
||||
|
||||
protected double[] chordLead = new double[DIVISIONS];
|
||||
protected double[] chordTrail = new double[DIVISIONS];
|
||||
protected double[] chordLength = new double[DIVISIONS];
|
||||
|
||||
private final double[] poly = new double[6];
|
||||
|
||||
private final double thickness;
|
||||
private final double bodyRadius;
|
||||
private final int finCount;
|
||||
private final double baseRotation;
|
||||
private final double cantAngle;
|
||||
protected final int interferenceFinCount;
|
||||
private final FinSet.CrossSection crossSection;
|
||||
|
||||
private final TubeFinSet tubeFinSet;
|
||||
|
||||
public TubeFinSetCalc(RocketComponent component) {
|
||||
super(component);
|
||||
if (!(component instanceof TubeFinSet)) {
|
||||
throw new IllegalArgumentException("Illegal component type " + component);
|
||||
}
|
||||
|
||||
tubeFinSet = (TubeFinSet) component;
|
||||
final TubeFinSet fin = (TubeFinSet) component;
|
||||
|
||||
thickness = fin.getThickness();
|
||||
bodyRadius = fin.getBodyRadius();
|
||||
finCount = 3 * fin.getFinCount();
|
||||
baseRotation = fin.getBaseRotation();
|
||||
cantAngle = 0;
|
||||
span = 2 * fin.getOuterRadius();
|
||||
finArea = span * fin.getLength();
|
||||
crossSection = FinSet.CrossSection.SQUARE;
|
||||
|
||||
calculateFinGeometry(fin);
|
||||
calculatePoly();
|
||||
interferenceFinCount = calculateInterferenceFinCount(fin);
|
||||
}
|
||||
|
||||
/*
|
||||
* Calculates the non-axial forces produced by the fins (normal and side forces,
|
||||
* pitch, yaw and roll moments, CP position, CNa).
|
||||
*/
|
||||
@Override
|
||||
public void calculateNonaxialForces(FlightConditions conditions, Transformation transform,
|
||||
AerodynamicForces forces) {
|
||||
|
||||
if (span < 0.001) {
|
||||
forces.setCm(0);
|
||||
forces.setCN(0);
|
||||
forces.setCNa(0);
|
||||
forces.setCP(Coordinate.NUL);
|
||||
forces.setCroll(0);
|
||||
forces.setCrollDamp(0);
|
||||
forces.setCrollForce(0);
|
||||
forces.setCside(0);
|
||||
forces.setCyaw(0);
|
||||
return;
|
||||
}
|
||||
|
||||
// Add warnings (radius/2 == diameter/4)
|
||||
if( (0 < bodyRadius) && (thickness > bodyRadius / 2)){
|
||||
warnings.add(Warning.THICK_FIN);
|
||||
}
|
||||
warnings.add(new Other("Tube fin support is experimental"));
|
||||
|
||||
//////// Calculate CNa. /////////
|
||||
|
||||
// One fin without interference (both sub- and supersonic):
|
||||
double cna1 = calculateFinCNa1(conditions);
|
||||
|
||||
// log.debug("Component cna1 = {}", cna1);
|
||||
|
||||
// Multiple fins with fin-fin interference
|
||||
double cna;
|
||||
double theta = conditions.getTheta();
|
||||
double angle = baseRotation + transform.getXrotation();
|
||||
|
||||
// Compute basic CNa without interference effects
|
||||
if (finCount == 1 || finCount == 2) {
|
||||
// Basic CNa from geometry
|
||||
double mul = 0;
|
||||
for (int i = 0; i < finCount; i++) {
|
||||
mul += MathUtil.pow2(Math.sin(theta - angle));
|
||||
angle += 2 * Math.PI / finCount;
|
||||
}
|
||||
cna = cna1 * mul;
|
||||
} else {
|
||||
// Basic CNa assuming full efficiency
|
||||
cna = cna1 * finCount / 2.0;
|
||||
}
|
||||
|
||||
// log.debug("Component cna = {}", cna);
|
||||
|
||||
// Take into account fin-fin interference effects
|
||||
switch (interferenceFinCount) {
|
||||
case 1:
|
||||
case 2:
|
||||
case 3:
|
||||
case 4:
|
||||
// No interference effect
|
||||
break;
|
||||
|
||||
case 5:
|
||||
cna *= 0.948;
|
||||
break;
|
||||
|
||||
case 6:
|
||||
cna *= 0.913;
|
||||
break;
|
||||
|
||||
case 7:
|
||||
cna *= 0.854;
|
||||
break;
|
||||
|
||||
case 8:
|
||||
cna *= 0.81;
|
||||
break;
|
||||
|
||||
default:
|
||||
// Assume 75% efficiency
|
||||
cna *= 0.75;
|
||||
break;
|
||||
}
|
||||
|
||||
|
||||
// Body-fin interference effect
|
||||
double r = bodyRadius;
|
||||
double tau = r / (span + r);
|
||||
if (Double.isNaN(tau) || Double.isInfinite(tau))
|
||||
tau = 0;
|
||||
cna *= 1 + tau; // Classical Barrowman
|
||||
// cna *= pow2(1 + tau); // Barrowman thesis (too optimistic??)
|
||||
// log.debug("Component cna = {}", cna);
|
||||
|
||||
// TODO: LOW: check for fin tip mach cone interference
|
||||
// (Barrowman thesis pdf-page 40)
|
||||
|
||||
// TODO: LOW: fin-fin mach cone effect, MIL-HDBK page 5-25
|
||||
|
||||
// Calculate CP position
|
||||
double x = macLead + calculateCPPos(conditions) * macLength;
|
||||
// log.debug("Component macLead = {}", macLead);
|
||||
// log.debug("Component macLength = {}", macLength);
|
||||
// log.debug("Component x = {}", x);
|
||||
|
||||
|
||||
// Calculate roll forces, reduce forcing above stall angle
|
||||
|
||||
// Without body-fin interference effect:
|
||||
// forces.CrollForce = fins * (macSpan+r) * cna1 * component.getCantAngle() /
|
||||
// conditions.getRefLength();
|
||||
// With body-fin interference effect:
|
||||
forces.setCrollForce(finCount * (macSpan + r) * cna1 * (1 + tau) * cantAngle / conditions.getRefLength());
|
||||
|
||||
if (conditions.getAOA() > STALL_ANGLE) {
|
||||
// System.out.println("Fin stalling in roll");
|
||||
forces.setCrollForce(forces.getCrollForce() * MathUtil.clamp(
|
||||
1 - (conditions.getAOA() - STALL_ANGLE) / (STALL_ANGLE / 2), 0, 1));
|
||||
}
|
||||
forces.setCrollDamp(calculateDampingMoment(conditions));
|
||||
forces.setCroll(forces.getCrollForce() - forces.getCrollDamp());
|
||||
|
||||
// System.out.printf(component.getName() + ": roll rate:%.3f force:%.3f damp:%.3f " +
|
||||
// "total:%.3f\n",
|
||||
// conditions.getRollRate(), forces.CrollForce, forces.CrollDamp, forces.Croll);
|
||||
|
||||
forces.setCNa(cna);
|
||||
forces.setCN(cna * MathUtil.min(conditions.getAOA(), STALL_ANGLE));
|
||||
forces.setCP(new Coordinate(x, 0, 0, cna));
|
||||
forces.setCm(forces.getCN() * x / conditions.getRefLength());
|
||||
|
||||
/*
|
||||
* TODO: HIGH: Compute actual side force and yaw moment.
|
||||
* This is not currently performed because it produces strange results for
|
||||
* stable rockets that have two fins in the front part of the fuselage,
|
||||
* where the rocket flies at an ever-increasing angle of attack. This may
|
||||
* be due to incorrect computation of pitch/yaw damping moments.
|
||||
*/
|
||||
// if (fins == 1 || fins == 2) {
|
||||
// forces.Cside = fins * cna1 * Math.cos(theta-angle) * Math.sin(theta-angle);
|
||||
// forces.Cyaw = fins * forces.Cside * x / conditions.getRefLength();
|
||||
// } else {
|
||||
// forces.Cside = 0;
|
||||
// forces.Cyaw = 0;
|
||||
// }
|
||||
forces.setCside(0);
|
||||
forces.setCyaw(0);
|
||||
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns the MAC length of the fin. This is required in the friction drag
|
||||
* computation.
|
||||
*
|
||||
* @return the MAC length of the fin.
|
||||
*/
|
||||
public double getMACLength() {
|
||||
return macLength;
|
||||
}
|
||||
|
||||
public double getMidchordPos() {
|
||||
return macLead + 0.5 * macLength;
|
||||
}
|
||||
|
||||
/**
|
||||
* Pre-calculates the fin geometry values.
|
||||
*/
|
||||
protected void calculateFinGeometry(TubeFinSet component) {
|
||||
|
||||
ar = 2 * pow2(span) / finArea;
|
||||
|
||||
Coordinate[] points = {
|
||||
Coordinate.NUL,
|
||||
new Coordinate(0, span),
|
||||
new Coordinate(component.getLength(), span),
|
||||
new Coordinate(component.getLength(), 0)
|
||||
};
|
||||
|
||||
|
||||
// Calculate the chord lead and trail positions and length
|
||||
|
||||
Arrays.fill(chordLead, Double.POSITIVE_INFINITY);
|
||||
Arrays.fill(chordTrail, Double.NEGATIVE_INFINITY);
|
||||
Arrays.fill(chordLength, 0);
|
||||
|
||||
for (int point = 1; point < points.length; point++) {
|
||||
double x1 = points[point - 1].x;
|
||||
double y1 = points[point - 1].y;
|
||||
double x2 = points[point].x;
|
||||
double y2 = points[point].y;
|
||||
|
||||
// Don't use the default EPSILON since it is too small
|
||||
// and causes too much numerical instability in the computation of x below
|
||||
if (MathUtil.equals(y1, y2, 0.001))
|
||||
continue;
|
||||
|
||||
int i1 = (int) (y1 * 1.0001 / span * (DIVISIONS - 1));
|
||||
int i2 = (int) (y2 * 1.0001 / span * (DIVISIONS - 1));
|
||||
i1 = MathUtil.clamp(i1, 0, DIVISIONS - 1);
|
||||
i2 = MathUtil.clamp(i2, 0, DIVISIONS - 1);
|
||||
if (i1 > i2) {
|
||||
int tmp = i2;
|
||||
i2 = i1;
|
||||
i1 = tmp;
|
||||
}
|
||||
|
||||
for (int i = i1; i <= i2; i++) {
|
||||
// Intersection point (x,y)
|
||||
double y = i * span / (DIVISIONS - 1);
|
||||
double x = (y - y2) / (y1 - y2) * x1 + (y1 - y) / (y1 - y2) * x2;
|
||||
if (x < chordLead[i])
|
||||
chordLead[i] = x;
|
||||
if (x > chordTrail[i])
|
||||
chordTrail[i] = x;
|
||||
|
||||
// TODO: LOW: If fin point exactly on chord line, might be counted twice:
|
||||
if (y1 < y2) {
|
||||
chordLength[i] -= x;
|
||||
} else {
|
||||
chordLength[i] += x;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Check and correct any inconsistencies
|
||||
for (int i = 0; i < DIVISIONS; i++) {
|
||||
if (Double.isInfinite(chordLead[i]) || Double.isInfinite(chordTrail[i]) ||
|
||||
Double.isNaN(chordLead[i]) || Double.isNaN(chordTrail[i])) {
|
||||
chordLead[i] = 0;
|
||||
chordTrail[i] = 0;
|
||||
}
|
||||
if (chordLength[i] < 0 || Double.isNaN(chordLength[i])) {
|
||||
chordLength[i] = 0;
|
||||
}
|
||||
if (chordLength[i] > chordTrail[i] - chordLead[i]) {
|
||||
chordLength[i] = chordTrail[i] - chordLead[i];
|
||||
}
|
||||
}
|
||||
|
||||
/* Calculate fin properties:
|
||||
*
|
||||
* macLength // MAC length
|
||||
* macLead // MAC leading edge position
|
||||
* macSpan // MAC spanwise position
|
||||
* ar // Fin aspect ratio (already set)
|
||||
* span // Fin span (already set)
|
||||
*/
|
||||
macLength = 0;
|
||||
macLead = 0;
|
||||
macSpan = 0;
|
||||
cosGamma = 0;
|
||||
cosGammaLead = 0;
|
||||
rollSum = 0;
|
||||
double area = 0;
|
||||
double radius = component.getBodyRadius();
|
||||
|
||||
final double dy = span / (DIVISIONS - 1);
|
||||
for (int i = 0; i < DIVISIONS; i++) {
|
||||
double length = chordTrail[i] - chordLead[i];
|
||||
double y = i * dy;
|
||||
|
||||
macLength += length * length;
|
||||
log.debug("macLength = {}, length = {}, i = {}", macLength, length, i);
|
||||
macSpan += y * length;
|
||||
macLead += chordLead[i] * length;
|
||||
area += length;
|
||||
rollSum += chordLength[i] * pow2(radius + y);
|
||||
|
||||
if (i > 0) {
|
||||
double dx = (chordTrail[i] + chordLead[i]) / 2 - (chordTrail[i - 1] + chordLead[i - 1]) / 2;
|
||||
cosGamma += dy / MathUtil.hypot(dx, dy);
|
||||
|
||||
dx = chordLead[i] - chordLead[i - 1];
|
||||
cosGammaLead += dy / MathUtil.hypot(dx, dy);
|
||||
}
|
||||
}
|
||||
|
||||
macLength *= dy;
|
||||
log.debug("macLength = {}", macLength);
|
||||
macSpan *= dy;
|
||||
macLead *= dy;
|
||||
area *= dy;
|
||||
rollSum *= dy;
|
||||
|
||||
macLength /= area;
|
||||
macSpan /= area;
|
||||
macLead /= area;
|
||||
cosGamma /= (DIVISIONS - 1);
|
||||
cosGammaLead /= (DIVISIONS - 1);
|
||||
}
|
||||
|
||||
/////////////// CNa1 calculation ////////////////
|
||||
|
||||
private static final double CNA_SUBSONIC = 0.9;
|
||||
private static final double CNA_SUPERSONIC = 1.5;
|
||||
private static final double CNA_SUPERSONIC_B = pow(pow2(CNA_SUPERSONIC) - 1, 1.5);
|
||||
private static final double GAMMA = 1.4;
|
||||
private static final LinearInterpolator K1, K2, K3;
|
||||
private static final PolyInterpolator cnaInterpolator = new PolyInterpolator(
|
||||
new double[] { CNA_SUBSONIC, CNA_SUPERSONIC },
|
||||
new double[] { CNA_SUBSONIC, CNA_SUPERSONIC },
|
||||
new double[] { CNA_SUBSONIC }
|
||||
);
|
||||
/* Pre-calculate the values for K1, K2 and K3 */
|
||||
static {
|
||||
// Up to Mach 5
|
||||
int n = (int) ((5.0 - CNA_SUPERSONIC) * 10);
|
||||
double[] x = new double[n];
|
||||
double[] k1 = new double[n];
|
||||
double[] k2 = new double[n];
|
||||
double[] k3 = new double[n];
|
||||
for (int i = 0; i < n; i++) {
|
||||
double M = CNA_SUPERSONIC + i * 0.1;
|
||||
double beta = MathUtil.safeSqrt(M * M - 1);
|
||||
x[i] = M;
|
||||
k1[i] = 2.0 / beta;
|
||||
k2[i] = ((GAMMA + 1) * pow(M, 4) - 4 * pow2(beta)) / (4 * pow(beta, 4));
|
||||
k3[i] = ((GAMMA + 1) * pow(M, 8) + (2 * pow2(GAMMA) - 7 * GAMMA - 5) * pow(M, 6) +
|
||||
10 * (GAMMA + 1) * pow(M, 4) + 8) / (6 * pow(beta, 7));
|
||||
}
|
||||
K1 = new LinearInterpolator(x, k1);
|
||||
K2 = new LinearInterpolator(x, k2);
|
||||
K3 = new LinearInterpolator(x, k3);
|
||||
|
||||
// System.out.println("K1[m="+CNA_SUPERSONIC+"] = "+k1[0]);
|
||||
// System.out.println("K2[m="+CNA_SUPERSONIC+"] = "+k2[0]);
|
||||
// System.out.println("K3[m="+CNA_SUPERSONIC+"] = "+k3[0]);
|
||||
}
|
||||
|
||||
protected double calculateFinCNa1(FlightConditions conditions) {
|
||||
double mach = conditions.getMach();
|
||||
double ref = conditions.getRefArea();
|
||||
double alpha = MathUtil.min(conditions.getAOA(),
|
||||
Math.PI - conditions.getAOA(), STALL_ANGLE);
|
||||
|
||||
// Subsonic case
|
||||
if (mach <= CNA_SUBSONIC) {
|
||||
return 2 * Math.PI * pow2(span) / (1 + MathUtil.safeSqrt(1 + (1 - pow2(mach)) *
|
||||
pow2(pow2(span) / (finArea * cosGamma)))) / ref;
|
||||
}
|
||||
|
||||
// Supersonic case
|
||||
if (mach >= CNA_SUPERSONIC) {
|
||||
return finArea * (K1.getValue(mach) + K2.getValue(mach) * alpha +
|
||||
K3.getValue(mach) * pow2(alpha)) / ref;
|
||||
}
|
||||
|
||||
// Transonic case, interpolate
|
||||
double subV, superV;
|
||||
double subD, superD;
|
||||
|
||||
double sq = MathUtil.safeSqrt(1 + (1 - pow2(CNA_SUBSONIC)) * pow2(span * span / (finArea * cosGamma)));
|
||||
subV = 2 * Math.PI * pow2(span) / ref / (1 + sq);
|
||||
subD = 2 * mach * Math.PI * pow(span, 6) / (pow2(finArea * cosGamma) * ref *
|
||||
sq * pow2(1 + sq));
|
||||
|
||||
superV = finArea * (K1.getValue(CNA_SUPERSONIC) + K2.getValue(CNA_SUPERSONIC) * alpha +
|
||||
K3.getValue(CNA_SUPERSONIC) * pow2(alpha)) / ref;
|
||||
superD = -finArea / ref * 2 * CNA_SUPERSONIC / CNA_SUPERSONIC_B;
|
||||
|
||||
// System.out.println("subV="+subV+" superV="+superV+" subD="+subD+" superD="+superD);
|
||||
|
||||
return cnaInterpolator.interpolate(mach, subV, superV, subD, superD, 0);
|
||||
}
|
||||
|
||||
private double calculateDampingMoment(FlightConditions conditions) {
|
||||
double rollRate = conditions.getRollRate();
|
||||
|
||||
if (Math.abs(rollRate) < 0.1)
|
||||
return 0;
|
||||
|
||||
double mach = conditions.getMach();
|
||||
double absRate = Math.abs(rollRate);
|
||||
|
||||
/*
|
||||
* At low speeds and relatively large roll rates (i.e. near apogee) the
|
||||
* fin tips rotate well above stall angle. In this case sum the chords
|
||||
* separately.
|
||||
*/
|
||||
if (absRate * (bodyRadius + span) / conditions.getVelocity() > 15 * Math.PI / 180) {
|
||||
double sum = 0;
|
||||
for (int i = 0; i < DIVISIONS; i++) {
|
||||
double dist = bodyRadius + span * i / DIVISIONS;
|
||||
double aoa = Math.min(absRate * dist / conditions.getVelocity(), 15 * Math.PI / 180);
|
||||
sum += chordLength[i] * dist * aoa;
|
||||
}
|
||||
sum = sum * (span / DIVISIONS) * 2 * Math.PI / conditions.getBeta() /
|
||||
(conditions.getRefArea() * conditions.getRefLength());
|
||||
|
||||
// System.out.println("SPECIAL: " +
|
||||
// (MathUtil.sign(rollRate) *component.getFinCount() * sum));
|
||||
return MathUtil.sign(rollRate) * finCount * sum;
|
||||
}
|
||||
|
||||
if (mach <= CNA_SUBSONIC) {
|
||||
// System.out.println("BASIC: "+
|
||||
// (component.getFinCount() * 2*Math.PI * rollRate * rollSum /
|
||||
// (conditions.getRefArea() * conditions.getRefLength() *
|
||||
// conditions.getVelocity() * conditions.getBeta())));
|
||||
|
||||
return finCount * 2 * Math.PI * rollRate * rollSum /
|
||||
(conditions.getRefArea() * conditions.getRefLength() *
|
||||
conditions.getVelocity() * conditions.getBeta());
|
||||
}
|
||||
if (mach >= CNA_SUPERSONIC) {
|
||||
|
||||
double vel = conditions.getVelocity();
|
||||
double k1 = K1.getValue(mach);
|
||||
double k2 = K2.getValue(mach);
|
||||
double k3 = K3.getValue(mach);
|
||||
|
||||
double sum = 0;
|
||||
|
||||
for (int i = 0; i < DIVISIONS; i++) {
|
||||
double y = i * span / (DIVISIONS - 1);
|
||||
double angle = rollRate * (bodyRadius + y) / vel;
|
||||
|
||||
sum += (k1 * angle + k2 * angle * angle + k3 * angle * angle * angle)
|
||||
* chordLength[i] * (bodyRadius + y);
|
||||
}
|
||||
|
||||
return finCount * sum * span / (DIVISIONS - 1) /
|
||||
(conditions.getRefArea() * conditions.getRefLength());
|
||||
}
|
||||
|
||||
// Transonic, do linear interpolation
|
||||
|
||||
FlightConditions cond = conditions.clone();
|
||||
cond.setMach(CNA_SUBSONIC - 0.01);
|
||||
double subsonic = calculateDampingMoment(cond);
|
||||
cond.setMach(CNA_SUPERSONIC + 0.01);
|
||||
double supersonic = calculateDampingMoment(cond);
|
||||
|
||||
return subsonic * (CNA_SUPERSONIC - mach) / (CNA_SUPERSONIC - CNA_SUBSONIC) +
|
||||
supersonic * (mach - CNA_SUBSONIC) / (CNA_SUPERSONIC - CNA_SUBSONIC);
|
||||
}
|
||||
|
||||
/**
|
||||
* Return the relative position of the CP along the mean aerodynamic chord.
|
||||
* Below mach 0.5 it is at the quarter chord, above mach 2 calculated using an
|
||||
* empirical formula, between these two using an interpolation polynomial.
|
||||
*
|
||||
* @param cond Mach speed used
|
||||
* @return CP position along the MAC
|
||||
*/
|
||||
private double calculateCPPos(FlightConditions cond) {
|
||||
double m = cond.getMach();
|
||||
// log.debug("m = {} ", m);
|
||||
if (m <= 0.5) {
|
||||
// At subsonic speeds CP at quarter chord
|
||||
return 0.25;
|
||||
}
|
||||
if (m >= 2) {
|
||||
// At supersonic speeds use empirical formula
|
||||
double beta = cond.getBeta();
|
||||
return (ar * beta - 0.67) / (2 * ar * beta - 1);
|
||||
}
|
||||
|
||||
// In between use interpolation polynomial
|
||||
double x = 1.0;
|
||||
double val = 0;
|
||||
|
||||
for (int i = 0; i < poly.length; i++) {
|
||||
val += poly[i] * x;
|
||||
x *= m;
|
||||
}
|
||||
// log.debug("val = {}", val);
|
||||
return val;
|
||||
}
|
||||
|
||||
/**
|
||||
* Calculate CP position interpolation polynomial coefficients from the
|
||||
* fin geometry. This is a fifth order polynomial that satisfies
|
||||
*
|
||||
* p(0.5)=0.25
|
||||
* p'(0.5)=0
|
||||
* p(2) = f(2)
|
||||
* p'(2) = f'(2)
|
||||
* p''(2) = 0
|
||||
* p'''(2) = 0
|
||||
*
|
||||
* where f(M) = (ar*sqrt(M^2-1) - 0.67) / (2*ar*sqrt(M^2-1) - 1).
|
||||
*
|
||||
* The values were calculated analytically in Mathematica. The coefficients
|
||||
* are used as poly[0] + poly[1]*x + poly[2]*x^2 + ...
|
||||
*/
|
||||
private void calculatePoly() {
|
||||
double denom = pow2(1 - 3.4641 * ar); // common denominator
|
||||
|
||||
poly[5] = (-1.58025 * (-0.728769 + ar) * (-0.192105 + ar)) / denom;
|
||||
poly[4] = (12.8395 * (-0.725688 + ar) * (-0.19292 + ar)) / denom;
|
||||
poly[3] = (-39.5062 * (-0.72074 + ar) * (-0.194245 + ar)) / denom;
|
||||
poly[2] = (55.3086 * (-0.711482 + ar) * (-0.196772 + ar)) / denom;
|
||||
poly[1] = (-31.6049 * (-0.705375 + ar) * (-0.198476 + ar)) / denom;
|
||||
poly[0] = (9.16049 * (-0.588838 + ar) * (-0.20624 + ar)) / denom;
|
||||
}
|
||||
|
||||
|
||||
// @SuppressWarnings("null")
|
||||
// public static void main(String arg[]) {
|
||||
// Rocket rocket = TestRocket.makeRocket();
|
||||
// FinSet finset = null;
|
||||
//
|
||||
// Iterator<RocketComponent> iter = rocket.deepIterator();
|
||||
// while (iter.hasNext()) {
|
||||
// RocketComponent c = iter.next();
|
||||
// if (c instanceof FinSet) {
|
||||
// finset = (FinSet)c;
|
||||
// break;
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// ((TrapezoidFinSet)finset).setHeight(0.10);
|
||||
// ((TrapezoidFinSet)finset).setRootChord(0.10);
|
||||
// ((TrapezoidFinSet)finset).setTipChord(0.10);
|
||||
// ((TrapezoidFinSet)finset).setSweep(0.0);
|
||||
//
|
||||
//
|
||||
// FinSetCalc calc = new FinSetCalc(finset);
|
||||
//
|
||||
// calc.calculateFinGeometry();
|
||||
// FlightConditions cond = new FlightConditions(new Configuration(rocket));
|
||||
// for (double m=0; m < 3; m+=0.05) {
|
||||
// cond.setMach(m);
|
||||
// cond.setAOA(0.0*Math.PI/180);
|
||||
// double cna = calc.calculateFinCNa1(cond);
|
||||
// System.out.printf("%5.2f "+cna+"\n", m);
|
||||
// }
|
||||
//
|
||||
// }
|
||||
|
||||
@Override
|
||||
public double calculateFrictionCD(FlightConditions conditions, double componentCf, WarningSet warnings) {
|
||||
if (tubeFinSet.getTubeSeparation() > MathUtil.EPSILON) {
|
||||
warnings.add(Warning.TUBE_SEPARATION);
|
||||
} else if (tubeFinSet.getTubeSeparation() < -MathUtil.EPSILON) {
|
||||
warnings.add(Warning.TUBE_OVERLAP);
|
||||
}
|
||||
|
||||
// Area of inner surface of tubes.
|
||||
final double innerArea = tubeFinSet.getFinCount() * tubeFinSet.getLength() * 2.0 * Math.PI * tubeFinSet.getInnerRadius();
|
||||
|
||||
// Area of the outer surface of tubes. Since roughly half
|
||||
// of the area is "masked" by the interstices between the tubes and the
|
||||
// body tube, only consider half the area
|
||||
final double outerArea = tubeFinSet.getFinCount() * tubeFinSet.getLength() * Math.PI * tubeFinSet.getOuterRadius();
|
||||
|
||||
// Area of the portion of the body tube masked by the tube fins
|
||||
final BodyTube parent = (BodyTube) tubeFinSet.getParent();
|
||||
final double maskedArea = tubeFinSet.getLength() * 2.0 * Math.PI * parent.getOuterRadius();
|
||||
|
||||
final double frictionCD = componentCf * (innerArea + outerArea - maskedArea);
|
||||
log.debug("frictionCD " + frictionCD);
|
||||
return frictionCD;
|
||||
}
|
||||
|
||||
@Override
|
||||
public double calculatePressureCD(FlightConditions conditions,
|
||||
double stagnationCD, double baseCD, WarningSet warnings) {
|
||||
|
||||
double mach = conditions.getMach();
|
||||
double cd = 0;
|
||||
|
||||
// Effective frontal area of tubes based on actual area and length of tubes
|
||||
|
||||
// Area of interstices. Assume they're long enough compared to other dimensions
|
||||
// they can be modeled as blunt
|
||||
|
||||
// Pressure fore-drag.
|
||||
if (crossSection == FinSet.CrossSection.AIRFOIL ||
|
||||
crossSection == FinSet.CrossSection.ROUNDED) {
|
||||
|
||||
// Round leading edge
|
||||
if (mach < 0.9) {
|
||||
cd = Math.pow(1 - pow2(mach), -0.417) - 1;
|
||||
} else if (mach < 1) {
|
||||
cd = 1 - 1.785 * (mach - 0.9);
|
||||
} else {
|
||||
cd = 1.214 - 0.502 / pow2(mach) + 0.1095 / pow2(pow2(mach));
|
||||
}
|
||||
|
||||
} else if (crossSection == FinSet.CrossSection.SQUARE) {
|
||||
cd = stagnationCD;
|
||||
} else {
|
||||
throw new UnsupportedOperationException("Unsupported fin profile: " + crossSection);
|
||||
}
|
||||
|
||||
// Slanted leading edge
|
||||
cd *= pow2(cosGammaLead);
|
||||
|
||||
// Trailing edge drag
|
||||
if (crossSection == FinSet.CrossSection.SQUARE) {
|
||||
cd += baseCD;
|
||||
} else if (crossSection == FinSet.CrossSection.ROUNDED) {
|
||||
cd += baseCD / 2;
|
||||
}
|
||||
// Airfoil assumed to have zero base drag
|
||||
|
||||
// Scale to correct reference area
|
||||
cd *= finCount * span * thickness / 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();
|
||||
}
|
||||
}
|
||||
11
core/src/net/sf/openrocket/rocketcomponent/Tube.java
Normal file
11
core/src/net/sf/openrocket/rocketcomponent/Tube.java
Normal file
@ -0,0 +1,11 @@
|
||||
package net.sf.openrocket.rocketcomponent;
|
||||
|
||||
import net.sf.openrocket.rocketcomponent.position.AxialMethod;
|
||||
|
||||
public abstract class Tube extends ExternalComponent implements Coaxial {
|
||||
|
||||
public Tube(AxialMethod relativePosition) {
|
||||
super(relativePosition);
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user