Made burn ctor not start burn.
[sw/motorsim] / 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         //Some constants to tune adaptive regression step\r
34         private static final double regStepIncreaseFactor = 1.01;\r
35         private static final double regStepDecreaseFactor = .5;\r
36         private static final Amount<Pressure> chamberPressureMaxDelta = Amount.valueOf(.5, SI.MEGA(SI.PASCAL));\r
37         \r
38         private static final Amount<Pressure> endPressure = Amount.valueOf(.1, RocketScience.PSI);\r
39         \r
40         private static Logger log = Logger.getLogger(Burn.class);\r
41         protected final Motor motor;\r
42         \r
43         private boolean burning = false;\r
44         private boolean done = false;\r
45         \r
46         public interface BurnProgressListener{\r
47                 public void setProgress(float p);\r
48         }\r
49         \r
50         private Set<BurnProgressListener> bpls = new HashSet<Burn.BurnProgressListener>();\r
51         \r
52         private static final Amount<Pressure> atmosphereicPressure = Amount.valueOf(101000, SI.PASCAL);\r
53         \r
54         public class Interval{\r
55                 public Amount<Duration> time;\r
56                 public Amount<Duration> dt;\r
57                 public Amount<Length> regression;\r
58                 public Amount<Pressure> chamberPressure;\r
59                 Amount<Mass> chamberProduct;\r
60                 public Amount<Force> thrust;\r
61 \r
62                 public String toString(){\r
63                         return time + " " + dt + " " + regression + " " + chamberPressure + " " + chamberProduct;\r
64                 }\r
65         }\r
66         \r
67         protected SortedMap<Amount<Duration>,Interval> data = new TreeMap<Amount<Duration>, Interval>();\r
68         \r
69         public SortedMap<Amount<Duration>,Interval> getData(){\r
70                 if ( !done )\r
71                         throw new IllegalStateException("Burn not complete!");\r
72                 return data;\r
73         }\r
74         \r
75         public Motor getMotor(){\r
76                 return motor;\r
77         }\r
78 \r
79         public Amount<Duration> burnTime(){\r
80                 return getData().lastKey();\r
81         }\r
82         \r
83         public Burn(Motor m){\r
84                 try {\r
85                         m.validate();\r
86                 } catch (ValidationException e) {\r
87                         throw new IllegalArgumentException("Invalid Motor: " + e.getMessage());\r
88                 }\r
89                 motor = m;\r
90         }\r
91         \r
92         public void addBurnProgressListener( BurnProgressListener bpl ){\r
93                 bpls.add(bpl);\r
94         }\r
95         \r
96         public void burn(){\r
97                 synchronized(this){\r
98                         if ( burning )\r
99                                 throw new IllegalStateException("Already burning!");\r
100                         burning = true;\r
101                 }\r
102                 log.info("Starting burn...");\r
103                 int endPressureSteps = 0;\r
104                 long start = new Date().getTime();\r
105                 \r
106                 Amount<Length> regStep = Amount.valueOf(0.01, SI.MILLIMETER);\r
107 \r
108                 Interval initial = new Interval();\r
109                 initial.time = Amount.valueOf(0, SI.SECOND);\r
110                 initial.dt = Amount.valueOf(0, SI.SECOND);\r
111                 initial.regression = Amount.valueOf(0, SI.MILLIMETER);\r
112                 initial.chamberPressure = atmosphereicPressure;\r
113                 initial.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
114                 initial.thrust = Amount.valueOf(0, SI.NEWTON);\r
115                 \r
116                 data.put(Amount.valueOf(0, SI.SECOND), initial);\r
117                 \r
118                 step:\r
119                 for ( int i = 0; i < 5000; i++ ) {\r
120                         assert(positive(regStep));\r
121                         regStep = regStep.times(regStepIncreaseFactor);\r
122                         \r
123                         Interval prev = data.get(data.lastKey());\r
124                         log.debug(prev);\r
125                         log.debug("Step " + i + " ==============================");\r
126                         Interval next = new Interval();\r
127                         \r
128                         Amount<Velocity> burnRate = motor.getFuel().burnRate(prev.chamberPressure);\r
129                         assert(positive(burnRate));\r
130                         \r
131                         log.debug("Burn Rate: " + burnRate);\r
132                         \r
133                         Amount<Duration> dt = regStep.divide(burnRate).to(Duration.UNIT);\r
134                         assert(positive(dt));\r
135                         next.dt = dt;\r
136                         \r
137 \r
138                         \r
139                         log.debug("Dt: " + dt);\r
140                         \r
141                         next.regression = prev.regression.plus(regStep);\r
142                         assert(positive(next.regression));\r
143                         \r
144                         log.debug("Regression: " + next.regression);\r
145                         \r
146                         //Update BurnProgressListeners\r
147                         Amount<Dimensionless> a = next.regression.divide(motor.getGrain().webThickness()).to(Dimensionless.UNIT);\r
148                         for (BurnProgressListener bpl : bpls ){\r
149                                 bpl.setProgress((float)a.doubleValue(Dimensionless.UNIT));\r
150                         }\r
151 \r
152                         \r
153                         next.time = prev.time.plus(dt);\r
154                         \r
155                         //log.debug("Vold: " + motor.getGrain().volume(prev.regression).to(SI.MILLIMETER.pow(3)));\r
156                         \r
157                         //log.debug("Vnew: " + motor.getGrain().volume(next.regression).to(SI.MILLIMETER.pow(3)));\r
158                         \r
159                         //TODO Amount<Volume> volumeBurnt = motor.getGrain().volume(prev.regression).minus(motor.getGrain().volume(next.regression));\r
160                         Amount<Volume> volumeBurnt = motor.getGrain().surfaceArea(prev.regression).times(regStep).to(Volume.UNIT);\r
161                         assert(positive(volumeBurnt));\r
162                         //log.info("Volume Burnt: " + volumeBurnt.to(SI.MILLIMETER.pow(3)));\r
163                         \r
164                         Amount<MassFlowRate> mGenRate = volumeBurnt.times(motor.getFuel().getIdealDensity().times(motor.getFuel().getDensityRatio())).divide(dt).to(MassFlowRate.UNIT);\r
165                         assert(positive(mGenRate));\r
166                         \r
167                         //log.debug("Mass Gen Rate: " + mGenRate);\r
168                         \r
169                         //Calculate specific gas constant\r
170                         Amount<?> specificGasConstant = Constants.R.divide(motor.getFuel().getCombustionProduct().getEffectiveMolarWeight());\r
171                         //This unit conversion helps JScience to convert nozzle flow rate to\r
172                         //kg/s a little later on I verified the conversion by hand and\r
173                         //JScience checks it too.\r
174                         specificGasConstant = convertSpecificGasConstantUnits(specificGasConstant);\r
175                         \r
176                         //Calculate chamber temperature\r
177                         Amount<Temperature> chamberTemp = motor.getFuel().getCombustionProduct().getIdealCombustionTemperature().times(motor.getFuel().getCombustionEfficiency());\r
178                         \r
179                         Amount<MassFlowRate> mNozzle;\r
180                         {\r
181                                 Amount<Pressure> pDiff = prev.chamberPressure.minus(atmosphereicPressure);\r
182                                 //log.debug("Pdiff: " + pDiff);\r
183                                 Amount<Area> aStar = motor.getNozzle().throatArea();\r
184                                 double k = motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats();\r
185                                 double kSide = Math.sqrt(k) * Math.pow((2/(k+1)) , (((k+1)/2)/(k-1)));\r
186                                 Amount<?> sqrtPart = specificGasConstant.times(chamberTemp).sqrt();\r
187                                 mNozzle = pDiff.times(aStar).times(kSide).divide(sqrtPart).to(MassFlowRate.UNIT);\r
188                                 //log.debug("Mass Exit Rate: " + mNozzle.to(MassFlowRate.UNIT));                \r
189                         }\r
190                         assert(positive(mNozzle));\r
191                         \r
192                         Amount<MassFlowRate> massStorageRate = mGenRate.minus(mNozzle);\r
193                         \r
194                         //log.debug("Mass Storage Rate: " + massStorageRate);\r
195 \r
196                         next.chamberProduct = prev.chamberProduct.plus(massStorageRate.times(dt));\r
197                         \r
198                         //Product can not go negative!\r
199                         if ( !positive(next.chamberProduct) ){\r
200                                 log.warn("ChamberProduct Negative on step " + i + "!, Adjusting regstep down and repeating step!");\r
201                                 regStep = regStep.times(regStepDecreaseFactor);\r
202                                 continue step;\r
203                         }\r
204                         assert(positive(next.chamberProduct));\r
205                         if ( next.chamberProduct.isLessThan(Amount.valueOf(0, SI.KILOGRAM)) )\r
206                                 next.chamberProduct = Amount.valueOf(0, SI.KILOGRAM);\r
207                         \r
208                         //log.debug("Chamber Product: " + next.chamberProduct);\r
209                         \r
210                         Amount<VolumetricDensity> combustionProductDensity = next.chamberProduct.divide(motor.getChamber().chamberVolume().minus(motor.getGrain().volume(next.regression))).to(VolumetricDensity.UNIT);\r
211                         \r
212                         log.debug("Product Density: " + combustionProductDensity);\r
213                         \r
214                         next.chamberPressure = combustionProductDensity.times(specificGasConstant).times(chamberTemp).plus(atmosphereicPressure).to(Pressure.UNIT);\r
215                         assert(positive(next.chamberPressure));\r
216                         \r
217                         next.chamberPressure = Amount.valueOf(\r
218                                         next.chamberPressure.doubleValue(SI.PASCAL),\r
219                                         SI.PASCAL);\r
220                         \r
221                         Amount<Pressure> dp = next.chamberPressure.minus(prev.chamberPressure);\r
222                         if ( dp.abs().isGreaterThan(chamberPressureMaxDelta)){\r
223                                 log.warn("DP " + dp + " too big!, Adjusting regstep down and repeating step!");\r
224                                 regStep = regStep.times(regStepDecreaseFactor);\r
225                                 continue step;\r
226                         }\r
227                         \r
228                         next.thrust = motor.getNozzle().thrust(next.chamberPressure, atmosphereicPressure, atmosphereicPressure, motor.getFuel().getCombustionProduct().getRatioOfSpecificHeats2Phase());\r
229                         assert(positive(next.thrust));\r
230                         \r
231                         if ( i > 100 && next.chamberPressure.minus(atmosphereicPressure).abs().isLessThan(endPressure)){\r
232                                 log.info("Pressure at ~Patm on step " + i);\r
233                                 endPressureSteps++;\r
234                                 if ( endPressureSteps > 5 )\r
235                                         break;\r
236                         }\r
237                         \r
238                         data.put(data.lastKey().plus(dt), next);\r
239                 }\r
240 \r
241                 long time = new Date().getTime() - start;\r
242                 log.info("Burn took " + time + " millis.");\r
243                 done = true;\r
244         }\r
245         \r
246         @SuppressWarnings("unchecked")\r
247         /*\r
248          * This converts the units of this constant to something JScience is able\r
249          * to work from. This conversion is unchecked at compile time, but\r
250          * JScience keeps me honest at runtime.\r
251          */\r
252         private Amount convertSpecificGasConstantUnits(Amount a){\r
253                 return a.to(\r
254                                 SI.METER.pow(2).divide(SI.SECOND.pow(2).times(SI.KELVIN)));\r
255         }\r
256         \r
257         public Amount<Pressure> pressure(Amount<Duration> time){\r
258                 return getData().get(time).chamberPressure;\r
259         }\r
260         \r
261         public Amount<Force> thrust(Amount<Duration> time){\r
262                 return getData().get(time).thrust;\r
263         }\r
264         \r
265         public Amount<Dimensionless> kn(Amount<Length> regression){\r
266                 return motor.getGrain().surfaceArea(regression).divide(motor.getNozzle().throatArea()).to(Dimensionless.UNIT);\r
267         }\r
268         \r
269         \r
270         private <Q extends Quantity> boolean positive(Amount<Q> a){\r
271                 return ( a.isGreaterThan(a.minus(a)) || a.equals(a.minus(a)));\r
272         }\r
273 \r
274 }\r