\r
\r
\r
+import java.util.Date;\r
+import java.util.HashSet;\r
+import java.util.Set;\r
import java.util.SortedMap;\r
import java.util.TreeMap;\r
\r
import javax.measure.quantity.Mass;\r
import javax.measure.quantity.MassFlowRate;\r
import javax.measure.quantity.Pressure;\r
+import javax.measure.quantity.Quantity;\r
import javax.measure.quantity.Temperature;\r
import javax.measure.quantity.Velocity;\r
import javax.measure.quantity.Volume;\r
import org.jscience.physics.amount.Amount;\r
import org.jscience.physics.amount.Constants;\r
\r
-import com.billkuker.rocketry.motorsim.fuel.KNSU;\r
-import com.billkuker.rocketry.motorsim.grain.BurnPanel;\r
-import com.billkuker.rocketry.motorsim.grain.CompoundGrain;\r
-import com.billkuker.rocketry.motorsim.grain.CoredCylindricalGrain;\r
-import com.billkuker.rocketry.motorsim.grain.ExtrudedGrain;\r
-import com.billkuker.rocketry.motorsim.grain.GrainPanel;\r
-import com.billkuker.rocketry.motorsim.grain.MultiGrain;\r
-import com.billkuker.rocketry.motorsim.visual.Chart;\r
+import com.billkuker.rocketry.motorsim.Validating.ValidationException;\r
\r
public class Burn {\r
- \r
private static Logger log = Logger.getLogger(Burn.class);\r
- protected final Motor motor;\r
\r
- private static final Amount<Pressure> atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL);\r
+ private static final BurnSettings settings = new BurnSettings();\r
+ public static final BurnSettings getBurnSettings(){\r
+ return settings;\r
+ }\r
\r
+ /**\r
+ * A class representing all the settigns one can change on a burn\r
+ * @author bkuker\r
+ */\r
+ public static class BurnSettings {\r
+ private BurnSettings(){};\r
+\r
+ public enum BurnVolumeMethod {\r
+ DeltaVolume,\r
+ SurfaceTimesRegression;\r
+ }\r
+ \r
+ private BurnVolumeMethod volumeMethod = BurnVolumeMethod.SurfaceTimesRegression;\r
+ private double regStepIncreaseFactor = 1.01;\r
+ private double regStepDecreaseFactor = .5;\r
+ private Amount<Pressure> chamberPressureMaxDelta = Amount.valueOf(.5, SI.MEGA(SI.PASCAL));\r
+ private Amount<Pressure> endPressure = Amount.valueOf(.1, RocketScience.PSI);\r
+ private Amount<Pressure> atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL);\r
+ \r
+ public void setVolumeMethod(BurnVolumeMethod volumeMethod) {\r
+ this.volumeMethod = volumeMethod;\r
+ }\r
+ public BurnVolumeMethod getVolumeMethod() {\r
+ return volumeMethod;\r
+ }\r
+ public double getRegStepIncreaseFactor() {\r
+ return regStepIncreaseFactor;\r
+ }\r
+ public void setRegStepIncreaseFactor(double regStepIncreaseFactor) {\r
+ this.regStepIncreaseFactor = regStepIncreaseFactor;\r
+ }\r
+ public double getRegStepDecreaseFactor() {\r
+ return regStepDecreaseFactor;\r
+ }\r
+ public void setRegStepDecreaseFactor(double regStepDecreaseFactor) {\r
+ this.regStepDecreaseFactor = regStepDecreaseFactor;\r
+ }\r
+ public Amount<Pressure> getChamberPressureMaxDelta() {\r
+ return chamberPressureMaxDelta;\r
+ }\r
+ public void setChamberPressureMaxDelta(Amount<Pressure> chamberPressureMaxDelta) {\r
+ this.chamberPressureMaxDelta = chamberPressureMaxDelta;\r
+ }\r
+ public Amount<Pressure> getEndPressure() {\r
+ return endPressure;\r
+ }\r
+ public void setEndPressure(Amount<Pressure> endPressure) {\r
+ this.endPressure = endPressure;\r
+ }\r
+ public Amount<Pressure> getAtmosphereicPressure() {\r
+ return atmosphereicPressure;\r
+ }\r
+ public void setAtmosphereicPressure(Amount<Pressure> atmosphereicPressure) {\r
+ this.atmosphereicPressure = atmosphereicPressure;\r
+ }\r
+ }\r
\r
- private static double combustionEfficency = 0.97;\r
+ protected final Motor motor;\r
\r
- private static double densityRatio = 0.96;\r
+ private boolean burning = false;\r
+ private boolean done = false;\r
\r
+ public interface BurnProgressListener{\r
+ public void setProgress(float p);\r
+ public void burnComplete();\r
+ }\r
+ \r
+ private Set<BurnProgressListener> bpls = new HashSet<Burn.BurnProgressListener>();\r
+\r
public class Interval{\r
- Amount<Duration> time;\r
+ public Amount<Duration> time;\r
+ public Amount<Duration> dt;\r
public Amount<Length> regression;\r
public Amount<Pressure> chamberPressure;\r
Amount<Mass> chamberProduct;\r
- Amount<Force> thrust;\r
+ public Amount<Force> thrust;\r
\r
public String toString(){\r
- return time + " " + regression + " " + chamberPressure + " " + chamberProduct;\r
+ return time + " " + dt + " " + regression + " " + chamberPressure + " " + chamberProduct;\r
}\r
}\r
\r
protected SortedMap<Amount<Duration>,Interval> data = new TreeMap<Amount<Duration>, Interval>();\r
\r
public SortedMap<Amount<Duration>,Interval> getData(){\r
+ if ( !done )\r
+ throw new IllegalStateException("Burn not complete!");\r
return data;\r
}\r
\r
}\r
\r
public Amount<Duration> burnTime(){\r
- return data.lastKey();\r
+ return getData().lastKey();\r
}\r
\r
public Burn(Motor m){\r
+ try {\r
+ m.validate();\r
+ } catch (ValidationException e) {\r
+ throw new IllegalArgumentException("Invalid Motor: " + e.getMessage());\r
+ }\r
motor = m;\r
}\r
\r
- private void burn(){\r
+ public void addBurnProgressListener( BurnProgressListener bpl ){\r
+ bpls.add(bpl);\r
+ if ( done )\r
+ bpl.burnComplete();\r
+ }\r
+ \r
+ public void burn(){\r
+ synchronized(this){\r
+ if ( burning )\r
+ throw new IllegalStateException("Already burning!");\r
+ burning = true;\r
+ }\r
log.info("Starting burn...");\r
+ for (BurnProgressListener bpl : bpls ){\r
+ bpl.setProgress(0);\r
+ }\r
\r
- Amount<Length> regStep = Amount.valueOf(0.0119904077, SI.MILLIMETER);\r
-\r
- //if ( motor.getGrain() instanceof Grain.DiscreteRegression )\r
- //regStep = ((Grain.DiscreteRegression)motor.getGrain()).optimalRegressionStep();\r
+ int endPressureSteps = 0;\r
+ long start = new Date().getTime();\r
\r
+ Amount<Length> regStep = Amount.valueOf(0.01, SI.MILLIMETER);\r
+\r
Interval initial = new Interval();\r
initial.time = Amount.valueOf(0, SI.SECOND);\r
+ initial.dt = Amount.valueOf(0, SI.SECOND);\r
initial.regression = Amount.valueOf(0, SI.MILLIMETER);\r
- initial.chamberPressure = atmosphereicPressure;\r
+ initial.chamberPressure = settings.getAtmosphereicPressure();\r
initial.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
initial.thrust = Amount.valueOf(0, SI.NEWTON);\r
\r
data.put(Amount.valueOf(0, SI.SECOND), initial);\r
\r
- for ( int i = 0; i < 5000; i++ ){\r
-\r
+ step:\r
+ for ( int i = 0; i < 5000; i++ ) {\r
+ assert(positive(regStep));\r
+ regStep = regStep.times(settings.getRegStepIncreaseFactor());\r
+ \r
Interval prev = data.get(data.lastKey());\r
log.debug(prev);\r
log.debug("Step " + i + " ==============================");\r
Interval next = new Interval();\r
\r
Amount<Velocity> burnRate = motor.getFuel().burnRate(prev.chamberPressure);\r
+ assert(positive(burnRate));\r
\r
log.debug("Burn Rate: " + burnRate);\r
\r
Amount<Duration> dt = regStep.divide(burnRate).to(Duration.UNIT);\r
+ assert(positive(dt));\r
+ next.dt = dt;\r
\r
- data.put(data.lastKey().plus(dt), next);\r
+\r
\r
log.debug("Dt: " + dt);\r
\r
next.regression = prev.regression.plus(regStep);\r
+ assert(positive(next.regression));\r
\r
- log.info("Regression: " + next.regression);\r
+ log.debug("Regression: " + next.regression);\r
+ \r
+ //Update BurnProgressListeners\r
+ Amount<Dimensionless> a = next.regression.divide(motor.getGrain().webThickness()).to(Dimensionless.UNIT);\r
+ for (BurnProgressListener bpl : bpls ){\r
+ bpl.setProgress((float)a.doubleValue(Dimensionless.UNIT));\r
+ }\r
+\r
\r
next.time = prev.time.plus(dt);\r
\r
- log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));\r
+ //log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));\r
\r
- log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));\r
+ //log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));\r
\r
- //TODO Amount<Volume> volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));\r
- Amount<Volume> volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);\r
+ Amount<Volume> volumeBurnt;\r
+ if ( settings.getVolumeMethod() == BurnSettings.BurnVolumeMethod.DeltaVolume ){\r
+ volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));\r
+ } else {\r
+ volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);\r
+ }\r
+ assert(positive(volumeBurnt));\r
+ //log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));\r
\r
- log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));\r
+ Amount<MassFlowRate> mGenRate = volumeBurnt.times(motor.getFuel().getIdealDensity().times(motor.getFuel().getDensityRatio())).divide(dt).to(MassFlowRate.UNIT);\r
+ assert(positive(mGenRate));\r
\r
- Amount<MassFlowRate> mGenRate = volumeBurnt.times(motor.getFuel().idealDensity().times(densityRatio)).divide(dt).to(MassFlowRate.UNIT);\r
+ //log.debug("Mass Gen Rate: " + mGenRate);\r
\r
- log.debug("Mass Gen Rate: " + mGenRate);\r
+ //Calculate specific gas constant\r
+ Amount<?> specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().getEffectiveMolarWeight());\r
+ //This unit conversion helps JScience to convert nozzle flow rate to\r
+ //kg/s a little later on I verified the conversion by hand and\r
+ //JScience checks it too.\r
+ specificGasConstant = convertSpecificGasConstantUnits(specificGasConstant);\r
\r
- Amount specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().effectiveMolarWeight());\r
- Amount<Temperature> chamberTemp = motor.getFuel().getCombustionProduct().idealCombustionTemperature().times(combustionEfficency);\r
+ //Calculate chamber temperature\r
+ Amount<Temperature> chamberTemp = motor.getFuel().getCombustionProduct().getIdealCombustionTemperature().times(motor.getFuel().getCombustionEfficiency());\r
\r
Amount<MassFlowRate> mNozzle;\r
{\r
- Amount<Pressure> pDiff = prev.chamberPressure.minus(atmosphereicPressure);\r
- \r
- //pDiff = Amount.valueOf(.7342, MPA).minus(atmosphereicPressure);\r
- \r
- log.debug("Pdiff: " + pDiff);\r
- \r
+ Amount<Pressure> pDiff = prev.chamberPressure.minus(settings.getAtmosphereicPressure());\r
+ //log.debug("Pdiff: " + pDiff);\r
Amount<Area> aStar = motor.getNozzle().throatArea();\r
- \r
- double k = motor.getFuel().getCombustionProduct().ratioOfSpecificHeats();\r
- \r
- log.debug("K: " + k);\r
- \r
- double kSide = Math.sqrt(k) * Math.pow((2/(k+1)) , (((k+1)/2)/(k-1))); //Math.pow(2/k+1, (k+1)/(2*(k-1)));\r
- \r
- log.debug("K-Part: (good)" + kSide);\r
- \r
- \r
- \r
- //This unit conversion helps JScience to convert nozzle flow rate to\r
- //kg/s a little later on I verified the conversion by hand and\r
- //JScience checks it too.\r
- specificGasConstant = specificGasConstant.to(\r
- SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN)));\r
- \r
- log.debug("Specific Gas Constant: (good)" + specificGasConstant);\r
- \r
- Amount sqrtPart = specificGasConstant.times(chamberTemp).sqrt();\r
-\r
- //Unit x = SI.JOULE.divide(SI.KILOGRAM).root(2);\r
- \r
- //sqrtPart = sqrtPart.times(Amount.valueOf(1, x));\r
- \r
- log.debug("Square Root Part: " + sqrtPart);\r
- \r
+ double k = motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats();\r
+ double kSide = Math.sqrt(k) * Math.pow((2/(k+1)) , (((k+1)/2)/(k-1)));\r
+ Amount<?> sqrtPart = specificGasConstant.times(chamberTemp).sqrt();\r
mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT);\r
- \r
- log.debug("Nozzle Flow: " + mNozzle);\r
- \r
- log.debug("Nozzle Flow: " + mNozzle.to(MassFlowRate.UNIT));\r
- \r
- \r
- \r
- \r
+ //log.debug("Mass Exit Rate: " + mNozzle.to(MassFlowRate.UNIT)); \r
}\r
+ assert(positive(mNozzle));\r
\r
Amount<MassFlowRate> massStorageRate = mGenRate.minus(mNozzle);\r
\r
- log.debug("Chamber Product rate: " + massStorageRate);\r
- \r
+ //log.debug("Mass Storage Rate: " + massStorageRate);\r
+\r
next.chamberProduct = prev.chamberProduct.plus(massStorageRate.times(dt));\r
\r
- log.debug("Chamber Product: " + next.chamberProduct);\r
+ //Product can not go negative!\r
+ if ( !positive(next.chamberProduct) ){\r
+ log.warn("ChamberProduct Negative on step " + i + "!, Adjusting regstep down and repeating step!");\r
+ regStep = regStep.times(settings.getRegStepDecreaseFactor());\r
+ continue step;\r
+ }\r
+ assert(positive(next.chamberProduct));\r
+ if ( next.chamberProduct.isLessThan(Amount.valueOf(0, SI.KILOGRAM)) )\r
+ next.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
+ \r
+ //log.debug("Chamber Product: " + next.chamberProduct);\r
\r
Amount<VolumetricDensity> combustionProductDensity = next.chamberProduct.divide(motor.getChamber().chamberVolume().minus(motor.getGrain().volume(next.regression))).to(VolumetricDensity.UNIT);\r
\r
log.debug("Product Density: " + combustionProductDensity);\r
\r
- next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(atmosphereicPressure).to(Pressure.UNIT);\r
+ next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(settings.getAtmosphereicPressure()).to(Pressure.UNIT);\r
+ assert(positive(next.chamberPressure));\r
\r
next.chamberPressure = Amount.valueOf(\r
next.chamberPressure.doubleValue(SI.PASCAL),\r
SI.PASCAL);\r
\r
- next.thrust = motor.getNozzle().thrust(next.chamberPressure, atmosphereicPressure, atmosphereicPressure, motor.getFuel().getCombustionProduct().ratioOfSpecificHeats2Phase());\r
+ Amount<Pressure> dp = next.chamberPressure.minus(prev.chamberPressure);\r
+ if ( dp.abs().isGreaterThan(settings.getChamberPressureMaxDelta())){\r
+ log.warn("DP " + dp + " too big!, Adjusting regstep down and repeating step!");\r
+ regStep = regStep.times(settings.getRegStepDecreaseFactor());\r
+ continue step;\r
+ }\r
+ \r
+ next.thrust = motor.getNozzle().thrust(next.chamberPressure, settings.getAtmosphereicPressure(), settings.getAtmosphereicPressure(), motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats2Phase());\r
+ assert(positive(next.thrust));\r
\r
- if ( next.chamberPressure.approximates(atmosphereicPressure)){\r
- log.info("Pressure at Patm on step " + i);\r
- break;\r
+ if ( i > 100 && next.chamberPressure.minus(settings.getAtmosphereicPressure()).abs().isLessThan(settings.getEndPressure())){\r
+ log.info("Pressure at ~Patm on step " + i);\r
+ endPressureSteps++;\r
+ if ( endPressureSteps > 5 )\r
+ break;\r
}\r
+ \r
+ data.put(data.lastKey().plus(dt), next);\r
}\r
- \r
+\r
+ long time = new Date().getTime() - start;\r
+ log.info("Burn took " + time + " millis.");\r
+ done = true;\r
+ for (BurnProgressListener bpl : bpls ){\r
+ bpl.burnComplete();\r
+ }\r
+ }\r
+ \r
+ @SuppressWarnings("unchecked")\r
+ /*\r
+ * This converts the units of this constant to something JScience is able\r
+ * to work from. This conversion is unchecked at compile time, but\r
+ * JScience keeps me honest at runtime.\r
+ */\r
+ private Amount convertSpecificGasConstantUnits(Amount a){\r
+ return a.to(\r
+ SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN)));\r
}\r
\r
public Amount<Pressure> pressure(Amount<Duration> time){\r
- return data.get(time).chamberPressure;\r
+ return getData().get(time).chamberPressure;\r
}\r
\r
public Amount<Force> thrust(Amount<Duration> time){\r
- return data.get(time).thrust;\r
+ return getData().get(time).thrust;\r
}\r
\r
public Amount<Dimensionless> kn(Amount<Length> regression){\r
return motor.getGrain().surfaceArea(regression).divide(motor.getNozzle().throatArea()).to(Dimensionless.UNIT);\r
}\r
\r
- public static void main( String args[]) throws Exception{\r
- Motor m = new Motor();\r
- m.setFuel(new KNSU());\r
- \r
- CylindricalChamber c = new CylindricalChamber();\r
- c.setLength(Amount.valueOf(200, SI.MILLIMETER));\r
- c.setID(Amount.valueOf(30, SI.MILLIMETER));\r
- m.setChamber(c);\r
- \r
- CoredCylindricalGrain g = new CoredCylindricalGrain();\r
- g.setLength(Amount.valueOf(70, SI.MILLIMETER));\r
- g.setOD(Amount.valueOf(30, SI.MILLIMETER));\r
- g.setID(Amount.valueOf(10, SI.MILLIMETER));\r
- m.setGrain(g);\r
- \r
- CoredCylindricalGrain g1 = new CoredCylindricalGrain();\r
- g1.setLength(Amount.valueOf(70, SI.MILLIMETER));\r
- g1.setOD(Amount.valueOf(30, SI.MILLIMETER));\r
- g1.setID(Amount.valueOf(18, SI.MILLIMETER));\r
- m.setGrain(g);\r
- \r
- CoredCylindricalGrain g2 = new CoredCylindricalGrain();\r
- g2.setLength(Amount.valueOf(70, SI.MILLIMETER));\r
- g2.setOD(Amount.valueOf(12, SI.MILLIMETER));\r
- g2.setID(Amount.valueOf(0, SI.MILLIMETER));\r
- g2.inhibit(true, false, false);\r
- m.setGrain(g);\r
- \r
- CompoundGrain cg = new CompoundGrain(g1, g2);\r
- \r
- m.setGrain( cg );\r
- \r
- //m.setGrain(new MultiGrain(g,2));\r
- \r
- //m.setGrain(new ExtrudedGrain());\r
- \r
- ConvergentDivergentNozzle n = new ConvergentDivergentNozzle();\r
- n.setThroatDiameter(Amount.valueOf(8.500, SI.MILLIMETER));\r
- n.setExitDiameter(Amount.valueOf(20.87, SI.MILLIMETER));\r
- n.setEfficiency(.87);\r
- m.setNozzle(n);\r
- \r
- Burn b = new Burn(m);\r
- \r
- b.burn();\r
- \r
- new BurnPanel(b).show();\r
- /*\r
- Chart<Duration, Pressure> r = new Chart<Duration, Pressure>(\r
- SI.SECOND,\r
- SI.MEGA(SI.PASCAL),\r
- b,\r
- "pressure");\r
- r.setDomain(b.data.keySet());\r
- r.show();\r
- \r
- Chart<Duration, Force> t = new Chart<Duration, Force>(\r
- SI.SECOND,\r
- SI.NEWTON,\r
- b,\r
- "thrust");\r
- t.setDomain(b.data.keySet());\r
- t.show();\r
- \r
- new GrainPanel( m.getGrain() ).show();*/\r
- \r
+ \r
+ private <Q extends Quantity> boolean positive(Amount<Q> a){\r
+ return ( a.isGreaterThan(a.minus(a)) || a.equals(a.minus(a)));\r
}\r
+\r
}\r
projects / sw / motorsim / blobdiff