git.gag.com Git - sw/motorsim/blob - src/com/billkuker/rocketry/motorsim/Burn.java

   1 package com.billkuker.rocketry.motorsim;\r
   2 \r
   3 \r
   4 \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
  10 \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
  25 \r
  26 import org.apache.log4j.Logger;\r
  27 import org.jscience.physics.amount.Amount;\r
  28 import org.jscience.physics.amount.Constants;\r
  29 \r
  30 import com.billkuker.rocketry.motorsim.Validating.ValidationException;\r
  31 \r
  32 public class Burn {\r
  33         private static Logger log = Logger.getLogger(Burn.class);\r
  34         \r
  35         private static final BurnSettings settings = new BurnSettings();\r
  36         public static final BurnSettings getBurnSettings(){\r
  37                 return settings;\r
  38         }\r
  39         \r
  40         /**\r
  41          * A class representing all the settigns one can change on a burn\r
  42          * @author bkuker\r
  43          */\r
  44         public static class BurnSettings {\r
  45                 private BurnSettings(){};\r
  46 \r
  47                 public enum BurnVolumeMethod {\r
  48                         DeltaVolume,\r
  49                         SurfaceTimesRegression;\r
  50                 }\r
  51                 \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
  58                 \r
  59                 public void setVolumeMethod(BurnVolumeMethod volumeMethod) {\r
  60                         this.volumeMethod = volumeMethod;\r
  61                 }\r
  62                 public BurnVolumeMethod getVolumeMethod() {\r
  63                         return volumeMethod;\r
  64                 }\r
  65                 public double getRegStepIncreaseFactor() {\r
  66                         return regStepIncreaseFactor;\r
  67                 }\r
  68                 public void setRegStepIncreaseFactor(double regStepIncreaseFactor) {\r
  69                         this.regStepIncreaseFactor = regStepIncreaseFactor;\r
  70                 }\r
  71                 public double getRegStepDecreaseFactor() {\r
  72                         return regStepDecreaseFactor;\r
  73                 }\r
  74                 public void setRegStepDecreaseFactor(double regStepDecreaseFactor) {\r
  75                         this.regStepDecreaseFactor = regStepDecreaseFactor;\r
  76                 }\r
  77                 public Amount<Pressure> getChamberPressureMaxDelta() {\r
  78                         return chamberPressureMaxDelta;\r
  79                 }\r
  80                 public void setChamberPressureMaxDelta(Amount<Pressure> chamberPressureMaxDelta) {\r
  81                         this.chamberPressureMaxDelta = chamberPressureMaxDelta;\r
  82                 }\r
  83                 public Amount<Pressure> getEndPressure() {\r
  84                         return endPressure;\r
  85                 }\r
  86                 public void setEndPressure(Amount<Pressure> endPressure) {\r
  87                         this.endPressure = endPressure;\r
  88                 }\r
  89                 public Amount<Pressure> getAtmosphereicPressure() {\r
  90                         return atmosphereicPressure;\r
  91                 }\r
  92                 public void setAtmosphereicPressure(Amount<Pressure> atmosphereicPressure) {\r
  93                         this.atmosphereicPressure = atmosphereicPressure;\r
  94                 }\r
  95         }\r
  96         \r
  97         protected final Motor motor;\r
  98         \r
  99         private boolean burning = false;\r
 100         private boolean done = false;\r
 101         \r
 102         public interface BurnProgressListener{\r
 103                 public void setProgress(float p);\r
 104                 public void burnComplete();\r
 105         }\r
 106                 \r
 107         private Set<BurnProgressListener> bpls = new HashSet<Burn.BurnProgressListener>();\r
 108 \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
 116                 public Amount<Mass> fuelBurnt;\r
 117 \r
 118                 public String toString(){\r
 119                         return time + " " + dt + " " + regression + " " + chamberPressure + " " + chamberProduct;\r
 120                 }\r
 121         }\r
 122         \r
 123         protected SortedMap<Amount<Duration>,Interval> data = new TreeMap<Amount<Duration>, Interval>();\r
 124         \r
 125         public SortedMap<Amount<Duration>,Interval> getData(){\r
 126                 if ( !done )\r
 127                         throw new IllegalStateException("Burn not complete!");\r
 128                 return data;\r
 129         }\r
 130         \r
 131         public Motor getMotor(){\r
 132                 return motor;\r
 133         }\r
 134 \r
 135         public Amount<Duration> burnTime(){\r
 136                 return getData().lastKey();\r
 137         }\r
 138         \r
 139         public Burn(Motor m){\r
 140                 try {\r
 141                         m.validate();\r
 142                 } catch (ValidationException e) {\r
 143                         throw new IllegalArgumentException("Invalid Motor: " + e.getMessage());\r
 144                 }\r
 145                 motor = m;\r
 146         }\r
 147         \r
 148         public void addBurnProgressListener( BurnProgressListener bpl ){\r
 149                 bpls.add(bpl);\r
 150                 if ( done )\r
 151                         bpl.burnComplete();\r
 152         }\r
 153         \r
 154         public void burn(){\r
 155                 synchronized(this){\r
 156                         if ( burning )\r
 157                                 throw new IllegalStateException("Already burning!");\r
 158                         burning = true;\r
 159                 }\r
 160                 log.info("Starting burn...");\r
 161                 for (BurnProgressListener bpl : bpls ){\r
 162                         bpl.setProgress(0);\r
 163                 }\r
 164                 \r
 165                 int endPressureSteps = 0;\r
 166                 long start = new Date().getTime();\r
 167                 \r
 168                 Amount<Length> regStep = Amount.valueOf(0.01, SI.MILLIMETER);\r
 169 \r
 170                 Interval initial = new Interval();\r
 171                 initial.time = Amount.valueOf(0, SI.SECOND);\r
 172                 initial.dt = Amount.valueOf(0, SI.SECOND);\r
 173                 initial.regression = Amount.valueOf(0, SI.MILLIMETER);\r
 174                 initial.chamberPressure = settings.getAtmosphereicPressure();\r
 175                 initial.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
 176                 initial.thrust = Amount.valueOf(0, SI.NEWTON);\r
 177                 initial.fuelBurnt = Amount.valueOf(0, SI.KILOGRAM);\r
 178                 \r
 179                 data.put(Amount.valueOf(0, SI.SECOND), initial);\r
 180                 \r
 181                 step:\r
 182                 for ( int i = 0; i < 5000; i++ ) {\r
 183                         assert(positive(regStep));\r
 184                         regStep = regStep.times(settings.getRegStepIncreaseFactor());\r
 185                         \r
 186                         Interval prev = data.get(data.lastKey());\r
 187                         log.debug(prev);\r
 188                         log.debug("Step " + i + " ==============================");\r
 189                         Interval next = new Interval();\r
 190                         \r
 191                         Amount<Velocity> burnRate = motor.getFuel().burnRate(prev.chamberPressure);\r
 192                         assert(positive(burnRate));\r
 193                         \r
 194                         log.debug("Burn Rate: " + burnRate);\r
 195                         \r
 196                         Amount<Duration> dt = regStep.divide(burnRate).to(Duration.UNIT);\r
 197                         assert(positive(dt));\r
 198                         next.dt = dt;\r
 199                         \r
 200 \r
 201                         \r
 202                         log.debug("Dt: " + dt);\r
 203                         \r
 204                         next.regression = prev.regression.plus(regStep);\r
 205                         assert(positive(next.regression));\r
 206                         \r
 207                         log.debug("Regression: " + next.regression);\r
 208                         \r
 209                         //Update BurnProgressListeners\r
 210                         Amount<Dimensionless> a = next.regression.divide(motor.getGrain().webThickness()).to(Dimensionless.UNIT);\r
 211                         for (BurnProgressListener bpl : bpls ){\r
 212                                 bpl.setProgress((float)a.doubleValue(Dimensionless.UNIT));\r
 213                         }\r
 214 \r
 215                         \r
 216                         next.time = prev.time.plus(dt);\r
 217                         \r
 218                         //log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));\r
 219                         \r
 220                         //log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));\r
 221                         \r
 222                         Amount<Volume> volumeBurnt;\r
 223                         if ( settings.getVolumeMethod() == BurnSettings.BurnVolumeMethod.DeltaVolume ){\r
 224                                 volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));\r
 225                         } else {\r
 226                                 volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);\r
 227                         }\r
 228                         assert(positive(volumeBurnt));\r
 229                         //log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));\r
 230                         \r
 231                         Amount<Mass> massBurnt = volumeBurnt.times(motor.getFuel().getIdealDensity().times(motor.getFuel().getDensityRatio())).to(Mass.UNIT);\r
 232                         next.fuelBurnt = massBurnt;\r
 233                         Amount<MassFlowRate> mGenRate = massBurnt.divide(dt).to(MassFlowRate.UNIT);\r
 234                         assert(positive(mGenRate));\r
 235                         \r
 236                         //log.debug("Mass Gen Rate: " + mGenRate);\r
 237                         \r
 238                         //Calculate specific gas constant\r
 239                         Amount<?> specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().getEffectiveMolarWeight());\r
 240                         //This unit conversion helps JScience to convert nozzle flow rate to\r
 241                         //kg/s a little later on I verified the conversion by hand and\r
 242                         //JScience checks it too.\r
 243                         specificGasConstant = convertSpecificGasConstantUnits(specificGasConstant);\r
 244                         \r
 245                         //Calculate chamber temperature\r
 246                         Amount<Temperature> chamberTemp = motor.getFuel().getCombustionProduct().getIdealCombustionTemperature().times(motor.getFuel().getCombustionEfficiency());\r
 247                         \r
 248                         Amount<MassFlowRate> mNozzle;\r
 249                         {\r
 250                                 Amount<Pressure> pDiff = prev.chamberPressure.minus(settings.getAtmosphereicPressure());\r
 251                                 //log.debug("Pdiff: " + pDiff);\r
 252                                 Amount<Area> aStar = motor.getNozzle().throatArea();\r
 253                                 double k = motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats();\r
 254                                 double kSide = Math.sqrt(k) * Math.pow((2/(k+1)) , (((k+1)/2)/(k-1)));\r
 255                                 Amount<?> sqrtPart = specificGasConstant.times(chamberTemp).sqrt();\r
 256                                 mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT);\r
 257                                 //log.debug("Mass Exit Rate: " + mNozzle.to(MassFlowRate.UNIT));                \r
 258                         }\r
 259                         assert(positive(mNozzle));\r
 260                         \r
 261                         Amount<MassFlowRate> massStorageRate = mGenRate.minus(mNozzle);\r
 262                         \r
 263                         //log.debug("Mass Storage Rate: " + massStorageRate);\r
 264 \r
 265                         next.chamberProduct = prev.chamberProduct.plus(massStorageRate.times(dt));\r
 266                         \r
 267                         //Product can not go negative!\r
 268                         if ( !positive(next.chamberProduct) ){\r
 269                                 log.warn("ChamberProduct Negative on step " + i + "!, Adjusting regstep down and repeating step!");\r
 270                                 regStep = regStep.times(settings.getRegStepDecreaseFactor());\r
 271                                 continue step;\r
 272                         }\r
 273                         assert(positive(next.chamberProduct));\r
 274                         if ( next.chamberProduct.isLessThan(Amount.valueOf(0, SI.KILOGRAM)) )\r
 275                                 next.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
 276                         \r
 277                         //log.debug("Chamber Product: " + next.chamberProduct);\r
 278                         \r
 279                         Amount<VolumetricDensity> combustionProductDensity = next.chamberProduct.divide(motor.getChamber().chamberVolume().minus(motor.getGrain().volume(next.regression))).to(VolumetricDensity.UNIT);\r
 280                         \r
 281                         log.debug("Product Density: " + combustionProductDensity);\r
 282                         \r
 283                         next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(settings.getAtmosphereicPressure()).to(Pressure.UNIT);\r
 284                         assert(positive(next.chamberPressure));\r
 285                         \r
 286                         next.chamberPressure = Amount.valueOf(\r
 287                                         next.chamberPressure.doubleValue(SI.PASCAL),\r
 288                                         SI.PASCAL);\r
 289                         \r
 290                         Amount<Pressure> dp = next.chamberPressure.minus(prev.chamberPressure);\r
 291                         if ( dp.abs().isGreaterThan(settings.getChamberPressureMaxDelta())){\r
 292                                 log.warn("DP " + dp + " too big!, Adjusting regstep down and repeating step!");\r
 293                                 regStep = regStep.times(settings.getRegStepDecreaseFactor());\r
 294                                 continue step;\r
 295                         }\r
 296                         \r
 297                         next.thrust = motor.getNozzle().thrust(next.chamberPressure, settings.getAtmosphereicPressure(), settings.getAtmosphereicPressure(), motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats2Phase());\r
 298                         assert(positive(next.thrust));\r
 299                         \r
 300                         if ( i > 100 && next.chamberPressure.minus(settings.getAtmosphereicPressure()).abs().isLessThan(settings.getEndPressure())){\r
 301                                 log.info("Pressure at ~Patm on step " + i);\r
 302                                 endPressureSteps++;\r
 303                                 if ( endPressureSteps > 5 )\r
 304                                         break;\r
 305                         }\r
 306                         \r
 307                         data.put(data.lastKey().plus(dt), next);\r
 308                 }\r
 309 \r
 310                 long time = new Date().getTime() - start;\r
 311                 log.info("Burn took " + time + " millis.");\r
 312                 done = true;\r
 313                 for (BurnProgressListener bpl : bpls ){\r
 314                         bpl.burnComplete();\r
 315                 }\r
 316         }\r
 317         \r
 318         @SuppressWarnings({ "unchecked", "rawtypes" })\r
 319         /*\r
 320          * This converts the units of this constant to something JScience is able\r
 321          * to work from. This conversion is unchecked at compile time, but\r
 322          * JScience keeps me honest at runtime.\r
 323          */\r
 324         private Amount convertSpecificGasConstantUnits(Amount a){\r
 325                 return a.to(\r
 326                                 SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN)));\r
 327         }\r
 328         \r
 329         public Amount<Pressure> pressure(Amount<Duration> time){\r
 330                 return getData().get(time).chamberPressure;\r
 331         }\r
 332         \r
 333         public Amount<Force> thrust(Amount<Duration> time){\r
 334                 return getData().get(time).thrust;\r
 335         }\r
 336         \r
 337         public Amount<Dimensionless> kn(Amount<Length> regression){\r
 338                 return motor.getGrain().surfaceArea(regression).divide(motor.getNozzle().throatArea()).to(Dimensionless.UNIT);\r
 339         }\r
 340         \r
 341         \r
 342         private <Q extends Quantity> boolean positive(Amount<Q> a){\r
 343                 return ( a.isGreaterThan(a.minus(a)) || a.equals(a.minus(a)));\r
 344         }\r
 345 \r
 346 }\r