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 import com.billkuker.rocketry.motorsim.fuel.KNSU;\r
  27 import com.billkuker.rocketry.motorsim.grain.BurnPanel;\r
  28 import com.billkuker.rocketry.motorsim.grain.CompoundGrain;\r
  29 import com.billkuker.rocketry.motorsim.grain.CoredCylindricalGrain;\r
  30 import com.billkuker.rocketry.motorsim.grain.ExtrudedGrain;\r
  31 import com.billkuker.rocketry.motorsim.grain.GrainPanel;\r
  32 import com.billkuker.rocketry.motorsim.grain.MultiGrain;\r
  33 import com.billkuker.rocketry.motorsim.visual.Chart;\r
  34 \r
  35 public class Burn {\r
  36         \r
  37         private static Logger log = Logger.getLogger(Burn.class);\r
  38         protected final Motor motor;\r
  39         \r
  40         private static final Amount<Pressure> atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL);\r
  41         \r
  42         \r
  43         private static double combustionEfficency = 0.97;\r
  44         \r
  45         private static double densityRatio = 0.96;\r
  46         \r
  47         public class Interval{\r
  48                 Amount<Duration> time;\r
  49                 public Amount<Length> regression;\r
  50                 public Amount<Pressure> chamberPressure;\r
  51                 Amount<Mass> chamberProduct;\r
  52                 Amount<Force> thrust;\r
  53 \r
  54                 public String toString(){\r
  55                         return time + " " + regression + " " + chamberPressure + " " + chamberProduct;\r
  56                 }\r
  57         }\r
  58         \r
  59         protected SortedMap<Amount<Duration>,Interval> data = new TreeMap<Amount<Duration>, Interval>();\r
  60         \r
  61         public SortedMap<Amount<Duration>,Interval> getData(){\r
  62                 return data;\r
  63         }\r
  64         \r
  65         public Motor getMotor(){\r
  66                 return motor;\r
  67         }\r
  68 \r
  69         public Amount<Duration> burnTime(){\r
  70                 return data.lastKey();\r
  71         }\r
  72         \r
  73         public Burn(Motor m){\r
  74                 motor = m;\r
  75         }\r
  76         \r
  77         private void burn(){\r
  78                 log.info("Starting burn...");\r
  79                 \r
  80                 Amount<Length> regStep = Amount.valueOf(0.0119904077, SI.MILLIMETER);\r
  81 \r
  82                 //if ( motor.getGrain() instanceof Grain.DiscreteRegression )\r
  83                         //regStep = ((Grain.DiscreteRegression)motor.getGrain()).optimalRegressionStep();\r
  84                 \r
  85                 Interval initial = new Interval();\r
  86                 initial.time = Amount.valueOf(0, SI.SECOND);\r
  87                 initial.regression = Amount.valueOf(0, SI.MILLIMETER);\r
  88                 initial.chamberPressure = atmosphereicPressure;\r
  89                 initial.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
  90                 initial.thrust = Amount.valueOf(0, SI.NEWTON);\r
  91                 \r
  92                 data.put(Amount.valueOf(0, SI.SECOND), initial);\r
  93                 \r
  94                 for ( int i = 0; i < 5000; i++ ){\r
  95 \r
  96                         Interval prev = data.get(data.lastKey());\r
  97                         log.debug(prev);\r
  98                         log.debug("Step " + i + " ==============================");\r
  99                         Interval next = new Interval();\r
 100                         \r
 101                         Amount<Velocity> burnRate = motor.getFuel().burnRate(prev.chamberPressure);\r
 102                         \r
 103                         log.debug("Burn Rate: " + burnRate);\r
 104                         \r
 105                         Amount<Duration> dt = regStep.divide(burnRate).to(Duration.UNIT);\r
 106                         \r
 107                         data.put(data.lastKey().plus(dt), next);\r
 108                         \r
 109                         log.debug("Dt: " + dt);\r
 110                         \r
 111                         next.regression = prev.regression.plus(regStep);\r
 112                         \r
 113                         log.info("Regression: " + next.regression);\r
 114                         \r
 115                         next.time = prev.time.plus(dt);\r
 116                         \r
 117                         log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));\r
 118                         \r
 119                         log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));\r
 120                         \r
 121                         //TODO Amount<Volume> volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));\r
 122                         Amount<Volume> volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);\r
 123                         \r
 124                         log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));\r
 125                         \r
 126                         Amount<MassFlowRate> mGenRate = volumeBurnt.times(motor.getFuel().idealDensity().times(densityRatio)).divide(dt).to(MassFlowRate.UNIT);\r
 127                         \r
 128                         log.debug("Mass Gen Rate: " + mGenRate);\r
 129                         \r
 130                         Amount specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().effectiveMolarWeight());\r
 131                         Amount<Temperature> chamberTemp = motor.getFuel().getCombustionProduct().idealCombustionTemperature().times(combustionEfficency);\r
 132                         \r
 133                         Amount<MassFlowRate> mNozzle;\r
 134                         {\r
 135                                 Amount<Pressure> pDiff = prev.chamberPressure.minus(atmosphereicPressure);\r
 136                                 \r
 137                                 //pDiff = Amount.valueOf(.7342, MPA).minus(atmosphereicPressure);\r
 138                                 \r
 139                                 log.debug("Pdiff: " + pDiff);\r
 140                                 \r
 141                                 Amount<Area> aStar = motor.getNozzle().throatArea();\r
 142                                 \r
 143                                 double k = motor.getFuel().getCombustionProduct().ratioOfSpecificHeats();\r
 144                                 \r
 145                                 log.debug("K: " + k);\r
 146                                 \r
 147                                 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
 148                                 \r
 149                                 log.debug("K-Part: (good)" + kSide);\r
 150                                 \r
 151                                 \r
 152                                 \r
 153                                 //This unit conversion helps JScience to convert nozzle flow rate to\r
 154                                 //kg/s a little later on I verified the conversion by hand and\r
 155                                 //JScience checks it too.\r
 156                                 specificGasConstant = specificGasConstant.to(\r
 157                                                 SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN)));\r
 158                                 \r
 159                                 log.debug("Specific Gas Constant: (good)" + specificGasConstant);\r
 160                                 \r
 161                                 Amount sqrtPart = specificGasConstant.times(chamberTemp).sqrt();\r
 162 \r
 163                                 //Unit x = SI.JOULE.divide(SI.KILOGRAM).root(2);\r
 164                                 \r
 165                                 //sqrtPart = sqrtPart.times(Amount.valueOf(1, x));\r
 166                                 \r
 167                                 log.debug("Square Root Part: " + sqrtPart);\r
 168                                 \r
 169                                 mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT);\r
 170                                 \r
 171                                 log.debug("Nozzle Flow: " + mNozzle);\r
 172                                 \r
 173                                 log.debug("Nozzle Flow: " + mNozzle.to(MassFlowRate.UNIT));\r
 174                                 \r
 175                                 \r
 176                                 \r
 177                                 \r
 178                         }\r
 179                         \r
 180                         Amount<MassFlowRate> massStorageRate = mGenRate.minus(mNozzle);\r
 181                         \r
 182                         log.debug("Chamber Product rate: " + massStorageRate);\r
 183                         \r
 184                         next.chamberProduct = prev.chamberProduct.plus(massStorageRate.times(dt));\r
 185                         \r
 186                         log.debug("Chamber Product: " + next.chamberProduct);\r
 187                         \r
 188                         Amount<VolumetricDensity> combustionProductDensity = next.chamberProduct.divide(motor.getChamber().chamberVolume().minus(motor.getGrain().volume(next.regression))).to(VolumetricDensity.UNIT);\r
 189                         \r
 190                         log.debug("Product Density: " + combustionProductDensity);\r
 191                         \r
 192                         next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(atmosphereicPressure).to(Pressure.UNIT);\r
 193                         \r
 194                         next.chamberPressure = Amount.valueOf(\r
 195                                         next.chamberPressure.doubleValue(SI.PASCAL),\r
 196                                         SI.PASCAL);\r
 197                         \r
 198                         next.thrust = motor.getNozzle().thrust(next.chamberPressure, atmosphereicPressure, atmosphereicPressure, motor.getFuel().getCombustionProduct().ratioOfSpecificHeats2Phase());\r
 199                         \r
 200                         if ( next.chamberPressure.approximates(atmosphereicPressure)){\r
 201                                 log.info("Pressure at Patm on step " + i);\r
 202                                 break;\r
 203                         }\r
 204                 }\r
 205                                 \r
 206         }\r
 207         \r
 208         public Amount<Pressure> pressure(Amount<Duration> time){\r
 209                 return data.get(time).chamberPressure;\r
 210         }\r
 211         \r
 212         public Amount<Force> thrust(Amount<Duration> time){\r
 213                 return data.get(time).thrust;\r
 214         }\r
 215         \r
 216         public Amount<Dimensionless> kn(Amount<Length> regression){\r
 217                 return motor.getGrain().surfaceArea(regression).divide(motor.getNozzle().throatArea()).to(Dimensionless.UNIT);\r
 218         }\r
 219         \r
 220         public static void main( String args[]) throws Exception{\r
 221                 Motor m = new Motor();\r
 222                 m.setFuel(new KNSU());\r
 223                 \r
 224                 CylindricalChamber c = new CylindricalChamber();\r
 225                 c.setLength(Amount.valueOf(200, SI.MILLIMETER));\r
 226                 c.setID(Amount.valueOf(30, SI.MILLIMETER));\r
 227                 m.setChamber(c);\r
 228                 \r
 229                 CoredCylindricalGrain g = new CoredCylindricalGrain();\r
 230                 g.setLength(Amount.valueOf(70, SI.MILLIMETER));\r
 231                 g.setOD(Amount.valueOf(30, SI.MILLIMETER));\r
 232                 g.setID(Amount.valueOf(10, SI.MILLIMETER));\r
 233                 m.setGrain(g);\r
 234                 \r
 235                 CoredCylindricalGrain g1 = new CoredCylindricalGrain();\r
 236                 g1.setLength(Amount.valueOf(70, SI.MILLIMETER));\r
 237                 g1.setOD(Amount.valueOf(30, SI.MILLIMETER));\r
 238                 g1.setID(Amount.valueOf(18, SI.MILLIMETER));\r
 239                 m.setGrain(g);\r
 240                 \r
 241                 CoredCylindricalGrain g2 = new CoredCylindricalGrain();\r
 242                 g2.setLength(Amount.valueOf(70, SI.MILLIMETER));\r
 243                 g2.setOD(Amount.valueOf(12, SI.MILLIMETER));\r
 244                 g2.setID(Amount.valueOf(0, SI.MILLIMETER));\r
 245                 g2.inhibit(true, false, false);\r
 246                 m.setGrain(g);\r
 247                 \r
 248                 CompoundGrain cg = new CompoundGrain(g1, g2);\r
 249                 \r
 250                 m.setGrain( cg );\r
 251                 \r
 252                 //m.setGrain(new MultiGrain(g,2));\r
 253                 \r
 254                 //m.setGrain(new ExtrudedGrain());\r
 255                 \r
 256                 ConvergentDivergentNozzle n = new ConvergentDivergentNozzle();\r
 257                 n.setThroatDiameter(Amount.valueOf(8.500, SI.MILLIMETER));\r
 258                 n.setExitDiameter(Amount.valueOf(20.87, SI.MILLIMETER));\r
 259                 n.setEfficiency(.87);\r
 260                 m.setNozzle(n);\r
 261                 \r
 262                 Burn b = new Burn(m);\r
 263                 \r
 264                 b.burn();\r
 265                 \r
 266                 new BurnPanel(b).show();\r
 267                 /*\r
 268                 Chart<Duration, Pressure> r = new Chart<Duration, Pressure>(\r
 269                                 SI.SECOND,\r
 270                                 SI.MEGA(SI.PASCAL),\r
 271                                 b,\r
 272                                 "pressure");\r
 273                 r.setDomain(b.data.keySet());\r
 274                 r.show();\r
 275                 \r
 276                 Chart<Duration, Force> t = new Chart<Duration, Force>(\r
 277                                 SI.SECOND,\r
 278                                 SI.NEWTON,\r
 279                                 b,\r
 280                                 "thrust");\r
 281                 t.setDomain(b.data.keySet());\r
 282                 t.show();\r
 283                 \r
 284                 new GrainPanel( m.getGrain() ).show();*/\r
 285                 \r
 286         }\r
 287 }\r