1 package com.billkuker.rocketry.motorsim;
\r
5 import java.util.Date;
\r
6 import java.util.HashSet;
\r
7 import java.util.Set;
\r
8 import java.util.SortedMap;
\r
9 import java.util.TreeMap;
\r
11 import javax.measure.quantity.Area;
\r
12 import javax.measure.quantity.Dimensionless;
\r
13 import javax.measure.quantity.Duration;
\r
14 import javax.measure.quantity.Force;
\r
15 import javax.measure.quantity.Length;
\r
16 import javax.measure.quantity.Mass;
\r
17 import javax.measure.quantity.MassFlowRate;
\r
18 import javax.measure.quantity.Pressure;
\r
19 import javax.measure.quantity.Quantity;
\r
20 import javax.measure.quantity.Temperature;
\r
21 import javax.measure.quantity.Velocity;
\r
22 import javax.measure.quantity.Volume;
\r
23 import javax.measure.quantity.VolumetricDensity;
\r
24 import javax.measure.unit.SI;
\r
26 import org.apache.log4j.Logger;
\r
27 import org.jscience.physics.amount.Amount;
\r
28 import org.jscience.physics.amount.Constants;
\r
30 import com.billkuker.rocketry.motorsim.Validating.ValidationException;
\r
33 private static Logger log = Logger.getLogger(Burn.class);
\r
35 private static final BurnSettings settings = new BurnSettings();
\r
36 public static final BurnSettings getBurnSettings(){
\r
41 * A class representing all the settigns one can change on a burn
\r
44 public static class BurnSettings {
\r
45 private BurnSettings(){};
\r
47 public enum BurnVolumeMethod {
\r
49 SurfaceTimesRegression;
\r
52 private BurnVolumeMethod volumeMethod = BurnVolumeMethod.SurfaceTimesRegression;
\r
53 private double regStepIncreaseFactor = 1.01;
\r
54 private double regStepDecreaseFactor = .5;
\r
55 private Amount<Pressure> chamberPressureMaxDelta = Amount.valueOf(.5, SI.MEGA(SI.PASCAL));
\r
56 private Amount<Pressure> endPressure = Amount.valueOf(.1, RocketScience.PSI);
\r
57 private Amount<Pressure> atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL);
\r
59 public void setVolumeMethod(BurnVolumeMethod volumeMethod) {
\r
60 this.volumeMethod = volumeMethod;
\r
62 public BurnVolumeMethod getVolumeMethod() {
\r
63 return volumeMethod;
\r
65 public double getRegStepIncreaseFactor() {
\r
66 return regStepIncreaseFactor;
\r
68 public void setRegStepIncreaseFactor(double regStepIncreaseFactor) {
\r
69 this.regStepIncreaseFactor = regStepIncreaseFactor;
\r
71 public double getRegStepDecreaseFactor() {
\r
72 return regStepDecreaseFactor;
\r
74 public void setRegStepDecreaseFactor(double regStepDecreaseFactor) {
\r
75 this.regStepDecreaseFactor = regStepDecreaseFactor;
\r
77 public Amount<Pressure> getChamberPressureMaxDelta() {
\r
78 return chamberPressureMaxDelta;
\r
80 public void setChamberPressureMaxDelta(Amount<Pressure> chamberPressureMaxDelta) {
\r
81 this.chamberPressureMaxDelta = chamberPressureMaxDelta;
\r
83 public Amount<Pressure> getEndPressure() {
\r
86 public void setEndPressure(Amount<Pressure> endPressure) {
\r
87 this.endPressure = endPressure;
\r
89 public Amount<Pressure> getAtmosphereicPressure() {
\r
90 return atmosphereicPressure;
\r
92 public void setAtmosphereicPressure(Amount<Pressure> atmosphereicPressure) {
\r
93 this.atmosphereicPressure = atmosphereicPressure;
\r
97 protected final Motor motor;
\r
99 private boolean burning = false;
\r
100 private boolean done = false;
\r
102 public interface BurnProgressListener{
\r
103 public void setProgress(float p);
\r
104 public void burnComplete();
\r
107 private Set<BurnProgressListener> bpls = new HashSet<Burn.BurnProgressListener>();
\r
109 public class Interval{
\r
110 public Amount<Duration> time;
\r
111 public Amount<Duration> dt;
\r
112 public Amount<Length> regression;
\r
113 public Amount<Pressure> chamberPressure;
\r
114 Amount<Mass> chamberProduct;
\r
115 public Amount<Force> thrust;
\r
117 public String toString(){
\r
118 return time + " " + dt + " " + regression + " " + chamberPressure + " " + chamberProduct;
\r
122 protected SortedMap<Amount<Duration>,Interval> data = new TreeMap<Amount<Duration>, Interval>();
\r
124 public SortedMap<Amount<Duration>,Interval> getData(){
\r
126 throw new IllegalStateException("Burn not complete!");
\r
130 public Motor getMotor(){
\r
134 public Amount<Duration> burnTime(){
\r
135 return getData().lastKey();
\r
138 public Burn(Motor m){
\r
141 } catch (ValidationException e) {
\r
142 throw new IllegalArgumentException("Invalid Motor: " + e.getMessage());
\r
147 public void addBurnProgressListener( BurnProgressListener bpl ){
\r
150 bpl.burnComplete();
\r
153 public void burn(){
\r
154 synchronized(this){
\r
156 throw new IllegalStateException("Already burning!");
\r
159 log.info("Starting burn...");
\r
160 for (BurnProgressListener bpl : bpls ){
\r
161 bpl.setProgress(0);
\r
164 int endPressureSteps = 0;
\r
165 long start = new Date().getTime();
\r
167 Amount<Length> regStep = Amount.valueOf(0.01, SI.MILLIMETER);
\r
169 Interval initial = new Interval();
\r
170 initial.time = Amount.valueOf(0, SI.SECOND);
\r
171 initial.dt = Amount.valueOf(0, SI.SECOND);
\r
172 initial.regression = Amount.valueOf(0, SI.MILLIMETER);
\r
173 initial.chamberPressure = settings.getAtmosphereicPressure();
\r
174 initial.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);
\r
175 initial.thrust = Amount.valueOf(0, SI.NEWTON);
\r
177 data.put(Amount.valueOf(0, SI.SECOND), initial);
\r
180 for ( int i = 0; i < 5000; i++ ) {
\r
181 assert(positive(regStep));
\r
182 regStep = regStep.times(settings.getRegStepIncreaseFactor());
\r
184 Interval prev = data.get(data.lastKey());
\r
186 log.debug("Step " + i + " ==============================");
\r
187 Interval next = new Interval();
\r
189 Amount<Velocity> burnRate = motor.getFuel().burnRate(prev.chamberPressure);
\r
190 assert(positive(burnRate));
\r
192 log.debug("Burn Rate: " + burnRate);
\r
194 Amount<Duration> dt = regStep.divide(burnRate).to(Duration.UNIT);
\r
195 assert(positive(dt));
\r
200 log.debug("Dt: " + dt);
\r
202 next.regression = prev.regression.plus(regStep);
\r
203 assert(positive(next.regression));
\r
205 log.debug("Regression: " + next.regression);
\r
207 //Update BurnProgressListeners
\r
208 Amount<Dimensionless> a = next.regression.divide(motor.getGrain().webThickness()).to(Dimensionless.UNIT);
\r
209 for (BurnProgressListener bpl : bpls ){
\r
210 bpl.setProgress((float)a.doubleValue(Dimensionless.UNIT));
\r
214 next.time = prev.time.plus(dt);
\r
216 //log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));
\r
218 //log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));
\r
220 Amount<Volume> volumeBurnt;
\r
221 if ( settings.getVolumeMethod() == BurnSettings.BurnVolumeMethod.DeltaVolume ){
\r
222 volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));
\r
224 volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);
\r
226 assert(positive(volumeBurnt));
\r
227 //log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));
\r
229 Amount<MassFlowRate> mGenRate = volumeBurnt.times(motor.getFuel().getIdealDensity().times(motor.getFuel().getDensityRatio())).divide(dt).to(MassFlowRate.UNIT);
\r
230 assert(positive(mGenRate));
\r
232 //log.debug("Mass Gen Rate: " + mGenRate);
\r
234 //Calculate specific gas constant
\r
235 Amount<?> specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().getEffectiveMolarWeight());
\r
236 //This unit conversion helps JScience to convert nozzle flow rate to
\r
237 //kg/s a little later on I verified the conversion by hand and
\r
238 //JScience checks it too.
\r
239 specificGasConstant = convertSpecificGasConstantUnits(specificGasConstant);
\r
241 //Calculate chamber temperature
\r
242 Amount<Temperature> chamberTemp = motor.getFuel().getCombustionProduct().getIdealCombustionTemperature().times(motor.getFuel().getCombustionEfficiency());
\r
244 Amount<MassFlowRate> mNozzle;
\r
246 Amount<Pressure> pDiff = prev.chamberPressure.minus(settings.getAtmosphereicPressure());
\r
247 //log.debug("Pdiff: " + pDiff);
\r
248 Amount<Area> aStar = motor.getNozzle().throatArea();
\r
249 double k = motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats();
\r
250 double kSide = Math.sqrt(k) * Math.pow((2/(k+1)) , (((k+1)/2)/(k-1)));
\r
251 Amount<?> sqrtPart = specificGasConstant.times(chamberTemp).sqrt();
\r
252 mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT);
\r
253 //log.debug("Mass Exit Rate: " + mNozzle.to(MassFlowRate.UNIT));
\r
255 assert(positive(mNozzle));
\r
257 Amount<MassFlowRate> massStorageRate = mGenRate.minus(mNozzle);
\r
259 //log.debug("Mass Storage Rate: " + massStorageRate);
\r
261 next.chamberProduct = prev.chamberProduct.plus(massStorageRate.times(dt));
\r
263 //Product can not go negative!
\r
264 if ( !positive(next.chamberProduct) ){
\r
265 log.warn("ChamberProduct Negative on step " + i + "!, Adjusting regstep down and repeating step!");
\r
266 regStep = regStep.times(settings.getRegStepDecreaseFactor());
\r
269 assert(positive(next.chamberProduct));
\r
270 if ( next.chamberProduct.isLessThan(Amount.valueOf(0, SI.KILOGRAM)) )
\r
271 next.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);
\r
273 //log.debug("Chamber Product: " + next.chamberProduct);
\r
275 Amount<VolumetricDensity> combustionProductDensity = next.chamberProduct.divide(motor.getChamber().chamberVolume().minus(motor.getGrain().volume(next.regression))).to(VolumetricDensity.UNIT);
\r
277 log.debug("Product Density: " + combustionProductDensity);
\r
279 next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(settings.getAtmosphereicPressure()).to(Pressure.UNIT);
\r
280 assert(positive(next.chamberPressure));
\r
282 next.chamberPressure = Amount.valueOf(
\r
283 next.chamberPressure.doubleValue(SI.PASCAL),
\r
286 Amount<Pressure> dp = next.chamberPressure.minus(prev.chamberPressure);
\r
287 if ( dp.abs().isGreaterThan(settings.getChamberPressureMaxDelta())){
\r
288 log.warn("DP " + dp + " too big!, Adjusting regstep down and repeating step!");
\r
289 regStep = regStep.times(settings.getRegStepDecreaseFactor());
\r
293 next.thrust = motor.getNozzle().thrust(next.chamberPressure, settings.getAtmosphereicPressure(), settings.getAtmosphereicPressure(), motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats2Phase());
\r
294 assert(positive(next.thrust));
\r
296 if ( i > 100 && next.chamberPressure.minus(settings.getAtmosphereicPressure()).abs().isLessThan(settings.getEndPressure())){
\r
297 log.info("Pressure at ~Patm on step " + i);
\r
298 endPressureSteps++;
\r
299 if ( endPressureSteps > 5 )
\r
303 data.put(data.lastKey().plus(dt), next);
\r
306 long time = new Date().getTime() - start;
\r
307 log.info("Burn took " + time + " millis.");
\r
309 for (BurnProgressListener bpl : bpls ){
\r
310 bpl.burnComplete();
\r
314 @SuppressWarnings("unchecked")
\r
316 * This converts the units of this constant to something JScience is able
\r
317 * to work from. This conversion is unchecked at compile time, but
\r
318 * JScience keeps me honest at runtime.
\r
320 private Amount convertSpecificGasConstantUnits(Amount a){
\r
322 SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN)));
\r
325 public Amount<Pressure> pressure(Amount<Duration> time){
\r
326 return getData().get(time).chamberPressure;
\r
329 public Amount<Force> thrust(Amount<Duration> time){
\r
330 return getData().get(time).thrust;
\r
333 public Amount<Dimensionless> kn(Amount<Length> regression){
\r
334 return motor.getGrain().surfaceArea(regression).divide(motor.getNozzle().throatArea()).to(Dimensionless.UNIT);
\r
338 private <Q extends Quantity> boolean positive(Amount<Q> a){
\r
339 return ( a.isGreaterThan(a.minus(a)) || a.equals(a.minus(a)));
\r