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.SortedMap;\r
   6 import java.util.TreeMap;\r
   7 \r
   8 import javax.measure.quantity.Area;\r
   9 import javax.measure.quantity.Dimensionless;\r
  10 import javax.measure.quantity.Duration;\r
  11 import javax.measure.quantity.Force;\r
  12 import javax.measure.quantity.Length;\r
  13 import javax.measure.quantity.Mass;\r
  14 import javax.measure.quantity.MassFlowRate;\r
  15 import javax.measure.quantity.Pressure;\r
  16 import javax.measure.quantity.Temperature;\r
  17 import javax.measure.quantity.Velocity;\r
  18 import javax.measure.quantity.Volume;\r
  19 import javax.measure.quantity.VolumetricDensity;\r
  20 import javax.measure.unit.SI;\r
  21 \r
  22 import org.apache.log4j.Logger;\r
  23 import org.jscience.physics.amount.Amount;\r
  24 import org.jscience.physics.amount.Constants;\r
  25 \r
  26 public class Burn {\r
  27         \r
  28         private static Logger log = Logger.getLogger(Burn.class);\r
  29         protected final Motor motor;\r
  30         \r
  31         private static final Amount<Pressure> atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL);\r
  32         \r
  33         public class Interval{\r
  34                 public Amount<Duration> time;\r
  35                 public Amount<Duration> dt;\r
  36                 public Amount<Length> regression;\r
  37                 public Amount<Pressure> chamberPressure;\r
  38                 Amount<Mass> chamberProduct;\r
  39                 public Amount<Force> thrust;\r
  40 \r
  41                 public String toString(){\r
  42                         return time + " " + dt + " " + regression + " " + chamberPressure + " " + chamberProduct;\r
  43                 }\r
  44         }\r
  45         \r
  46         protected SortedMap<Amount<Duration>,Interval> data = new TreeMap<Amount<Duration>, Interval>();\r
  47         \r
  48         public SortedMap<Amount<Duration>,Interval> getData(){\r
  49                 return data;\r
  50         }\r
  51         \r
  52         public Motor getMotor(){\r
  53                 return motor;\r
  54         }\r
  55 \r
  56         public Amount<Duration> burnTime(){\r
  57                 return data.lastKey();\r
  58         }\r
  59         \r
  60         public Burn(Motor m){\r
  61                 motor = m;\r
  62                 burn();\r
  63         }\r
  64         \r
  65         private void burn(){\r
  66                 log.info("Starting burn...");\r
  67                 \r
  68                 Amount<Length> regStep = Amount.valueOf(0.0119904077, SI.MILLIMETER);\r
  69 \r
  70                 //if ( motor.getGrain() instanceof Grain.DiscreteRegression )\r
  71                         //regStep = ((Grain.DiscreteRegression)motor.getGrain()).optimalRegressionStep();\r
  72                 \r
  73                 Interval initial = new Interval();\r
  74                 initial.time = Amount.valueOf(0, SI.SECOND);\r
  75                 initial.dt = Amount.valueOf(0, SI.SECOND);\r
  76                 initial.regression = Amount.valueOf(0, SI.MILLIMETER);\r
  77                 initial.chamberPressure = atmosphereicPressure;\r
  78                 initial.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
  79                 initial.thrust = Amount.valueOf(0, SI.NEWTON);\r
  80                 \r
  81                 data.put(Amount.valueOf(0, SI.SECOND), initial);\r
  82                 \r
  83                 for ( int i = 0; i < 5000; i++ ){\r
  84 \r
  85                         Interval prev = data.get(data.lastKey());\r
  86                         log.debug(prev);\r
  87                         log.debug("Step " + i + " ==============================");\r
  88                         Interval next = new Interval();\r
  89                         \r
  90                         Amount<Velocity> burnRate = motor.getFuel().burnRate(prev.chamberPressure);\r
  91                         \r
  92                         log.debug("Burn Rate: " + burnRate);\r
  93                         \r
  94                         Amount<Duration> dt = regStep.divide(burnRate).to(Duration.UNIT);\r
  95                         next.dt = dt;\r
  96                         \r
  97                         data.put(data.lastKey().plus(dt), next);\r
  98                         \r
  99                         log.debug("Dt: " + dt);\r
 100                         \r
 101                         next.regression = prev.regression.plus(regStep);\r
 102                         \r
 103                         log.info("Regression: " + next.regression);\r
 104                         \r
 105                         next.time = prev.time.plus(dt);\r
 106                         \r
 107                         log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));\r
 108                         \r
 109                         log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));\r
 110                         \r
 111                         //TODO Amount<Volume> volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));\r
 112                         Amount<Volume> volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);\r
 113                         log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));\r
 114                         \r
 115                         Amount<MassFlowRate> mGenRate = volumeBurnt.times(motor.getFuel().getIdealDensity().times(motor.getFuel().getDensityRatio())).divide(dt).to(MassFlowRate.UNIT);\r
 116                         log.debug("Mass Gen Rate: " + mGenRate);\r
 117                         \r
 118                         //Calculate specific gas constant\r
 119                         Amount specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().getEffectiveMolarWeight());\r
 120                         //This unit conversion helps JScience to convert nozzle flow rate to\r
 121                         //kg/s a little later on I verified the conversion by hand and\r
 122                         //JScience checks it too.\r
 123                         specificGasConstant = convertSpecificGasConstantUnits(specificGasConstant);\r
 124                         \r
 125                         //Calculate chamber temperature\r
 126                         Amount<Temperature> chamberTemp = motor.getFuel().getCombustionProduct().getIdealCombustionTemperature().times(motor.getFuel().getCombustionEfficiency());\r
 127                         \r
 128                         Amount<MassFlowRate> mNozzle;\r
 129                         {\r
 130                                 Amount<Pressure> pDiff = prev.chamberPressure.minus(atmosphereicPressure);\r
 131                                 log.debug("Pdiff: " + pDiff);\r
 132                                 Amount<Area> aStar = motor.getNozzle().throatArea();\r
 133                                 double k = motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats();\r
 134                                 double kSide = Math.sqrt(k) * Math.pow((2/(k+1)) , (((k+1)/2)/(k-1)));\r
 135                                 Amount sqrtPart = specificGasConstant.times(chamberTemp).sqrt();\r
 136                                 mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT);\r
 137                                 log.debug("Mass Exit Rate: " + mNozzle.to(MassFlowRate.UNIT));          \r
 138                         }\r
 139                         \r
 140                         Amount<MassFlowRate> massStorageRate = mGenRate.minus(mNozzle);\r
 141                         \r
 142                         log.debug("Mass Storage Rate: " + massStorageRate);\r
 143 \r
 144                         next.chamberProduct = prev.chamberProduct.plus(massStorageRate.times(dt));\r
 145                         \r
 146                         //Product can not go negative!\r
 147                         if ( next.chamberProduct.isLessThan(Amount.valueOf(0, SI.KILOGRAM)) )\r
 148                                 next.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
 149                         \r
 150                         log.debug("Chamber Product: " + next.chamberProduct);\r
 151                         \r
 152                         Amount<VolumetricDensity> combustionProductDensity = next.chamberProduct.divide(motor.getChamber().chamberVolume().minus(motor.getGrain().volume(next.regression))).to(VolumetricDensity.UNIT);\r
 153                         \r
 154                         log.debug("Product Density: " + combustionProductDensity);\r
 155                         \r
 156                         next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(atmosphereicPressure).to(Pressure.UNIT);\r
 157                         \r
 158                         next.chamberPressure = Amount.valueOf(\r
 159                                         next.chamberPressure.doubleValue(SI.PASCAL),\r
 160                                         SI.PASCAL);\r
 161                         \r
 162                         next.thrust = motor.getNozzle().thrust(next.chamberPressure, atmosphereicPressure, atmosphereicPressure, motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats2Phase());\r
 163                         \r
 164                         if ( i > 100 && next.chamberPressure.approximates(atmosphereicPressure)){\r
 165                                 log.info("Pressure at Patm on step " + i);\r
 166                                 break;\r
 167                         }\r
 168                 }\r
 169 \r
 170         }\r
 171         \r
 172         @SuppressWarnings("unchecked")\r
 173         /*\r
 174          * This converts the units of this constant to something JScience is able\r
 175          * to work from. This conversion is unchecked at compile time, but\r
 176          * JScience keeps me honest at runtime.\r
 177          */\r
 178         private Amount convertSpecificGasConstantUnits(Amount a){\r
 179                 return a.to(\r
 180                                 SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN)));\r
 181         }\r
 182         \r
 183         public Amount<Pressure> pressure(Amount<Duration> time){\r
 184                 return data.get(time).chamberPressure;\r
 185         }\r
 186         \r
 187         public Amount<Force> thrust(Amount<Duration> time){\r
 188                 return data.get(time).thrust;\r
 189         }\r
 190         \r
 191         public Amount<Dimensionless> kn(Amount<Length> regression){\r
 192                 return motor.getGrain().surfaceArea(regression).divide(motor.getNozzle().throatArea()).to(Dimensionless.UNIT);\r
 193         }\r
 194         \r
 195 }\r