hongshaorou/openrocket
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Modified integrate() to use geometry of actual frustums in geometric calculations instead of cylinder approximations. This turned into a pretty complete rewrite of integrate(), and the creation of several helper functions.

Browse Source 
Changed the name of integrate() to calculateProperties() to better express the purpose rather than the method of calculation. It's looking forward to a future PR when I'll be modifying BodyTube to use SymmetricComponent's variable caching.
This commit is contained in:
JoePfeiffer 2023-12-07 09:17:08 -07:00
parent e92e2bc3a9
commit 3e1341271c
1 changed files with 208 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

312
core/src/net/sf/openrocket/rocketcomponent/SymmetricComponent.java
View File

@ -1,7 +1,5 @@
package net.sf.openrocket.rocketcomponent;
import static net.sf.openrocket.util.MathUtil.pow2;
import java.util.ArrayList;
import java.util.Collection;
import java.util.List;
@ -11,6 +9,8 @@ import net.sf.openrocket.rocketcomponent.position.AxialMethod;
import net.sf.openrocket.util.BoundingBox;
import net.sf.openrocket.util.Coordinate;
import net.sf.openrocket.util.MathUtil;
import static net.sf.openrocket.util.MathUtil.pow2;
/**
* Class for an axially symmetric rocket component generated by rotating
@ -27,7 +27,6 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
protected boolean filled = false;
protected double thickness = DEFAULT_THICKNESS;
// Cached data, default values signify not calculated
private double wetArea = Double.NaN;
@ -38,7 +37,6 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
private double longitudinalInertia = Double.NaN;
private double rotationalInertia = Double.NaN;
private Coordinate cg = null;
public SymmetricComponent() {
super();
@ -171,7 +169,6 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
clearPreset();
}
/**
* Adds component bounds at a number of points between 0...length.
*/
@ -185,8 +182,6 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
}
return list;
}
@Override
protected void loadFromPreset(ComponentPreset preset) {
@ -200,94 +195,73 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
super.loadFromPreset(preset);
}
/**
* Calculate volume of the component by integrating over the length of the component.
* The method caches the result, so subsequent calls are instant. Subclasses may
* override this method for simple shapes and use this method as necessary.
* Obtain volume of component
*
* @return The volume of the component.
*/
@Override
public double getComponentVolume() {
if (Double.isNaN(volume)) {
integrate();
calculateProperties();
}
return volume;
}
/**
* Calculate full (filled) volume of the component by integrating over the length
* of the component. The method caches the result, so subsequent calls are instant.
* Subclasses may override this method for simple shapes and use this method as
* necessary.
* Obtain full (filled) volume of the component
*
* @return The filled volume of the component.
*/
public double getFullVolume() {
if (Double.isNaN(fullVolume)) {
integrate();
calculateProperties();
}
return fullVolume;
}
/**
* Calculate the wetted area of the component by integrating over the length
* of the component. The method caches the result, so subsequent calls are instant.
* Subclasses may override this method for simple shapes and use this method as
* necessary.
*
* Obtain wetted area of the component
*
* @return The wetted area of the component.
*/
public double getComponentWetArea() {
if (Double.isNaN(wetArea)) {
integrate();
calculateProperties();
}
return wetArea;
}
/**
* Calculate the planform area of the component by integrating over the length of
* the component. The method caches the result, so subsequent calls are instant.
* Subclasses may override this method for simple shapes and use this method as
* necessary.
* Obtain planform area of the component
*
* @return The planform area of the component.
*/
public double getComponentPlanformArea() {
if (Double.isNaN(planArea)) {
integrate();
calculateProperties();
}
return planArea;
}
/**
* Calculate the planform center X-coordinate of the component by integrating over
* the length of the component. The planform center is defined as
* Obtain planform area of the component, defined as
* <pre> integrate(x*2*r(x)) / planform area </pre>
* The method caches the result, so subsequent calls are instant. Subclasses may
* override this method for simple shapes and use this method as necessary.
*
* @return The planform center of the component.
*/
public double getComponentPlanformCenter() {
if (Double.isNaN(planCenter)) {
integrate();
calculateProperties();
}
return planCenter;
}
/**
* Return the CG and mass of the component. Subclasses may
* Obtain the CG and mass of the component. Subclasses may
* override this method for simple shapes and use this method as necessary.
* Note that subclasses will typically include shoulders in their calculations
*
* @return The CG+mass of the component.
*/
@ -297,47 +271,158 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
}
/**
* Calculate CG of the symmetric component by integrating over the length of the component.
* The method caches the result, so subsequent calls are instant. We need this method because subclasses
* override getComponentCG() and include mass of shoulders
* Obtain the CG and mass of the symmetric component. We need this method because subclasses
* override getComponentCG() and include mass of shoulders (if any), while moment calculations in
* this class assume no shoulders
*
* @return The CG+mass of the component.
* @return The CG+mass of the symmetric component.
*/
private Coordinate getSymmetricComponentCG() {
if (cg == null)
integrate();
calculateProperties();
return cg;
}
/*
* Obtain longitudinal unit inertia Ixx
*
* @return longitudinal unit inertia Ixx
*/
@Override
public double getLongitudinalUnitInertia() {
if (Double.isNaN(longitudinalInertia)) {
integrate();
calculateProperties();
}
return longitudinalInertia;
}
/*
* Obtain rotational unit inertia Iyy = Izz
*
* @return rotational unit inertia
*/
@Override
public double getRotationalUnitInertia() {
if (Double.isNaN(rotationalInertia)) {
integrate();
calculateProperties();
}
return rotationalInertia;
}
/**
* helper method for calculateProperties()
* returns a Coordinate with the CG (relative to fore end of frustum) and volume of a filled conical
* frustum. Result is also correct for cases of r1=r2, r1=0, and r2=0
* Caution! actually returns 3/PI times correct value, avoiding extra operations in loop. Gets corrected at end
* of numerical integration loop
* @param l length (height) of frustum
* @param r1 radius of fore end of frustum
* @param r2 radius of aft end of frustum
* @return Coordinate containing volume (as mass) and CG of frustum
*/
private Coordinate calculateCG(double l, double r1, double r2) {
final double volume = l * (pow2(r1) + r1 * r2 + pow2(r2));
final double cg;
if (volume < MathUtil.EPSILON) {
cg = l/2.0;
} else {
cg = l * (pow2(r1) + 2.0 * r1 * r2 + 3*pow2(r2)) / (4.0 * (pow2(r1) + r1*r2 + pow2(r2)));
}
return new Coordinate(cg, 0, 0, volume);
}
/**
* helper method for calculateProperties()
* returns the rotational moment of inertia of a solid frustum
* handles special case of cylinder correctly
* @param l length (height) of frustum
* @param r1 radius of fore end of frustum
* @param r2 radius of aft end of frustum
* @return rotational moment of inertia
*/
private double calculateRotMOI(double l, double r1, double r2) {
// check for cylinder special case
if (Math.abs(r1 - r2) < MathUtil.EPSILON) {
return pow2(r1)/2.0;
}
return 3.0 * (Math.pow(r2, 5) - Math.pow(r1, 5)) / (10.0 * (Math.pow(r2, 3) - Math.pow(r1, 3)));
}
/**
* helper method for calculateLongMOI(). Computes longitudinal MOI of a cone
* uses 'h' for height of cone to avoid confusion with x1 and x2 in calculateProperties()
* requires multiplication by pi later
* @param h height of cone
* @param r radius of cone
* @return longitudinal moment of inertia relative to tip of cone
*/
private double calculateLongMOICone(double h, double r) {
final double m = pow2(r) * h;
final double Ixx = 3*m*(pow2(r)/20.0 + pow2(h)/5.0);
return Ixx;
}
/**
* helper method for calculateProperties(). Computes longitudinal MOI of a solid frustum relative to CG
* calculates by subtracting MOI of a shorter cone from MOI of longer cone
* more readable than deriving and using the final formula
* handles special case of cylinder
* Caution! Requires multiplication by pi later
* @param l length of frustum
* @param r1 radius of forend of frustum
* @param r2 radius of aft end of frustum
* @param cg Coordinate with CG and mass of frustum (relative to fore end of frustum)
* @return longitudinal MOI
*/
private double calculateLongMOI(double l, double r1, double r2, Coordinate cg) {
// check for cylinder special case
double moi;
if (Math.abs(r1 - r2) < MathUtil.EPSILON) {
// compute MOI of cylinder relative to CG of cylinder
moi = cg.weight * (3 * pow2(r1) + pow2(l))/12.0;
return moi;
}
// is the frustum "small end forward" or "small end aft"?
double shiftCG = cg.x;
if (r1 > r2) {
final double tmp = r1;
r1 = r2;
r2 = tmp;
shiftCG = l - cg.x;
}
// Find the heights of the two cones. Note that the h1 and h2 being calculated here
// are NOT the x1 and x2 used in calculateProperties()
final double h2 = l*r2/(r2 - r1);
final double h1 = h2*r1/r2;
final double moi1 = calculateLongMOICone(h1, r1);
final double moi2 = calculateLongMOICone(h2, r2);
// compute MOI relative to tip of cones (they share the same tip, of course)
moi = moi2 - moi1;
// use parallel axis theorem to move MOI to be relative to CG of frustum.
moi = moi - pow2(h1 + shiftCG) * cg.weight;
return moi;
}
/**
* Performs integration over the length of the component and updates the cached variables.
*/
private void integrate() {
private void calculateProperties() {
wetArea = 0;
planArea = 0;
planCenter = 0;
fullVolume = 0;
volume = 0;
cg = Coordinate.NUL;
longitudinalInertia = 0;
rotationalInertia = 0;
@ -351,68 +436,89 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
// Integrate for volume, CG, wetted area, planform area, and moments of inertia
for (int n = 0; n < DIVISIONS; n++) {
/*
* x1 and x2 are the bounds on this section
* x1 and x2 are the bounds on this division
* hyp is the length of the hypotenuse from r1 to r2
* height is the y-axis height of the component if not filled
* r1o and r2o are the outer radii
* r1i and r2i are the inner radii
*/
// get x bounds and length for this division
final double x1 = n * getLength() / DIVISIONS;
final double r1 = getRadius(x1);
final double x2 = (n + 1) * getLength() / DIVISIONS;
final double r2 = getRadius(x2);
final double l = x2 - x1;
final double xmean = (x1 + x2) / 2;
final double outer = (r1 + r2) / 2;
final double hyp = MathUtil.hypot(r2 - r1, l);
// get outer and inner radii
final double r1o = getRadius(x1);
final double r2o = getRadius(x2);
// use thickness and angle of outer wall to get height of ring
final double hyp = MathUtil.hypot(r2o - r1o, l);
final double height = thickness * hyp / l;
// Volume differential elements
final double inner;
final double dV;
final double dFullV;
dFullV = Math.PI * l * (r1 * r1 + r1 * r2 + r2 * r2) / 3.0;
if ( filled || r1 < height || r2 < height ) {
inner = 0;
dV = dFullV;
// get inner radii.
final double r1i;
final double r2i;
if (filled) {
r1i = 0;
r2i = 0;
} else {
// hollow
// Thickness is normal to the surface of the component
// here we use simple trig to project the Thickness
// on to the y dimension (radius).
if (r1 < height || r2 < height) {
// Filled portion of piece
dV = dFullV;
inner = 0;
} else {
// Hollow portion of piece
dV = MathUtil.max(Math.PI* l * height * (r1 + r2 - height), 0);
inner = Math.max(outer - height, 0.);
}
// Tiny inaccuracy is introduced on a division where one end is closed and other is open.
r1i = MathUtil.max(r1o - height, 0);
r2i = MathUtil.max(r2o - height, 0);
}
// find volume and CG of (possibly hollow) frustum
final Coordinate fullCG = calculateCG(l, r1o, r2o);
final Coordinate innerCG = calculateCG(l, r1i, r2i);
final double dFullV = fullCG.weight;
final double dV = fullCG.weight - innerCG.weight;
final double dCG = (fullCG.x*fullCG.weight - innerCG.x*innerCG.weight)/dV;
// First moment, used later for CG calculation
final double dCGx = dV * (x1 + dCG);
// rotational moment of inertia
final double Izzo = calculateRotMOI(l, r1o, r2o);
final double Izzi = calculateRotMOI(l, r1i, r2i);
final double Izz = Izzo * fullCG.weight - Izzi * innerCG.weight;
// longitudinal moment of inertia -- axis through CG of division
double Ixx = calculateLongMOI(l, r1o, r2o, fullCG) - calculateLongMOI(l, r1i, r2i, innerCG);
// move to axis through forward end of component
Ixx += dV * pow2(x1 + dCG);
// Add to the volume-related components
volume += dV;
fullVolume += dFullV;
cgx += xmean * dV;
// Wetted area ( * PI at the end)
wetArea += hyp * (r1 + r2);
// Planform area & center
final double p = l * (r1 + r2);
planArea += p;
planCenter += xmean * p;
cgx += dCGx;
rotationalInertia += Izz;
longitudinalInertia += Ixx;
rotationalInertia += dV * (pow2(outer) + pow2(inner)) / 2;
longitudinalInertia += dV * ((3 * (pow2(outer) + pow2(inner)) + pow2(l)) / 12 + pow2(xmean));
// Wetted area ( * PI at the end)
wetArea += (r1o + r2o) * Math.sqrt(pow2(r1o - r2o) + pow2(l));
// Planform area & moment
final double dA = l*(r1o + r2o);
planArea += dA;
final double planMoment = dA*x1 + 2.0*pow2(l)*(r1o/6.0 + r2o/3.0);
planCenter += planMoment;
}
wetArea *= Math.PI;
// get unit moments of inertia
rotationalInertia /= volume;
longitudinalInertia /= volume;
if (planArea > 0)
planCenter /= planArea;
// Correct for deferred constant factors
volume *= Math.PI / 3.0;
fullVolume *= Math.PI / 3.0;
cgx *= Math.PI / 3.0;
wetArea *= Math.PI;
if (volume < 0.0000000001) { // 0.1 mm^3
volume = 0;
@ -424,18 +530,16 @@ public abstract class SymmetricComponent extends BodyComponent implements BoxBou
cg = new Coordinate(cgx / volume, 0, 0, getMaterial().getDensity() * volume );
}
// a part so small it has no volume can't contribute to moment of inertia
// a component so small it has no volume can't contribute to moment of inertia
if (MathUtil.equals(volume, 0)) {
rotationalInertia = 0;
longitudinalInertia = 0;
return;
}
rotationalInertia /= volume;
longitudinalInertia /= volume;
// Shift longitudinal inertia to CG
longitudinalInertia = longitudinalInertia - pow2(getSymmetricComponentCG().x);
longitudinalInertia = longitudinalInertia - pow2(cg.x);
}
/**